AppFemPy.cpp

Go to the documentation of this file.
00001 /***************************************************************************
00002  *   Copyright (c) 2008 Jürgen Riegel (juergen.riegel@web.de)              *
00003  *                                                                         *
00004  *   This file is part of the FreeCAD CAx development system.              *
00005  *                                                                         *
00006  *   This library is free software; you can redistribute it and/or         *
00007  *   modify it under the terms of the GNU Library General Public           *
00008  *   License as published by the Free Software Foundation; either          *
00009  *   version 2 of the License, or (at your option) any later version.      *
00010  *                                                                         *
00011  *   This library  is distributed in the hope that it will be useful,      *
00012  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00013  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00014  *   GNU Library General Public License for more details.                  *
00015  *                                                                         *
00016  *   You should have received a copy of the GNU Library General Public     *
00017  *   License along with this library; see the file COPYING.LIB. If not,    *
00018  *   write to the Free Software Foundation, Inc., 59 Temple Place,         *
00019  *   Suite 330, Boston, MA  02111-1307, USA                                *
00020  *                                                                         *
00021  ***************************************************************************/
00022 
00023 
00024 #include "PreCompiled.h"
00025 #ifndef _PreComp_
00026 # include <Python.h>
00027 # include <memory>
00028 #endif
00029 
00030 #include <Base/Console.h>
00031 #include <Base/VectorPy.h>
00032 #include <Base/PlacementPy.h>
00033 #include <App/Application.h>
00034 #include <App/Document.h>
00035 #include <App/DocumentObject.h>
00036 #include <App/DocumentObjectPy.h>
00037 #include <Mod/Mesh/App/Core/MeshKernel.h>
00038 #include <Mod/Mesh/App/Core/Evaluation.h>
00039 #include <Mod/Mesh/App/Core/Iterator.h>
00040 
00041 #include <SMESH_Gen.hxx>
00042 #include <SMESH_Group.hxx>
00043 #include <SMESHDS_Mesh.hxx>
00044 #include <SMDS_MeshNode.hxx>
00045 #include <StdMeshers_MaxLength.hxx>
00046 #include <StdMeshers_LocalLength.hxx>
00047 #include <StdMeshers_NumberOfSegments.hxx>
00048 #include <StdMeshers_AutomaticLength.hxx>
00049 #include <StdMeshers_TrianglePreference.hxx>
00050 #include <StdMeshers_MEFISTO_2D.hxx>
00051 #include <StdMeshers_Deflection1D.hxx>
00052 #include <StdMeshers_MaxElementArea.hxx>
00053 #include <StdMeshers_Regular_1D.hxx>
00054 #include <StdMeshers_QuadranglePreference.hxx>
00055 #include <StdMeshers_Quadrangle_2D.hxx>
00056 
00057 #include <StdMeshers_LengthFromEdges.hxx>
00058 #include <StdMeshers_NotConformAllowed.hxx>
00059 #include <StdMeshers_Arithmetic1D.hxx>
00060 
00061 #include "FemMesh.h"
00062 #include "FemMeshObject.h"
00063 #include "FemMeshPy.h"
00064 
00065 #include <cstdlib>
00066 
00067 #include <Standard_Real.hxx>
00068 #include "Base/Vector3D.h"
00069 
00070 using namespace Fem;
00071 
00072 
00073 /* module functions */
00074 static PyObject * read(PyObject *self, PyObject *args)
00075 {
00076     const char* Name;
00077     if (!PyArg_ParseTuple(args, "s",&Name))
00078         return NULL;
00079 
00080     PY_TRY {
00081         std::auto_ptr<FemMesh> mesh(new FemMesh);
00082         try {
00083             mesh->read(Name);
00084             return new FemMeshPy(mesh.release());
00085         }
00086         catch(...) {
00087             PyErr_SetString(PyExc_Exception, "Loading of mesh was aborted");
00088             return NULL;
00089         }
00090     } PY_CATCH;
00091 
00092     Py_Return;
00093 }
00094 
00095 static PyObject * open(PyObject *self, PyObject *args)
00096 {
00097     const char* Name;
00098     if (!PyArg_ParseTuple(args, "s",&Name))
00099         return NULL;
00100 
00101     PY_TRY {
00102         std::auto_ptr<FemMesh> mesh(new FemMesh);
00103         mesh->read(Name);
00104         Base::FileInfo file(Name);
00105         // create new document and add Import feature
00106         App::Document *pcDoc = App::GetApplication().newDocument("Unnamed");
00107         FemMeshObject *pcFeature = static_cast<FemMeshObject *>
00108             (pcDoc->addObject("Fem::FemMeshObject", file.fileNamePure().c_str()));
00109         pcFeature->Label.setValue(file.fileNamePure().c_str());
00110         pcFeature->FemMesh.setValuePtr(mesh.get());
00111         (void)mesh.release();
00112         pcFeature->purgeTouched();
00113     } PY_CATCH;
00114 
00115     Py_Return;
00116 }
00117 
00118 
00119 static PyObject * SMESH_PCA(PyObject *self, PyObject *args)
00120 {
00121         PyObject *input;
00122 
00123         if (!PyArg_ParseTuple(args, "O",&input))
00124                 return NULL;
00125 
00126         PY_TRY {
00127 
00128                 FemMeshPy *inputMesh = static_cast<FemMeshPy*>(input); 
00129                 MeshCore::MeshKernel aMesh;
00130                 MeshCore::MeshPointArray vertices;
00131                 vertices.clear();
00132                 MeshCore::MeshFacetArray faces;
00133                 faces.clear();
00134                 MeshCore::MeshPoint current_node;
00135                 SMDS_NodeIteratorPtr aNodeIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->nodesIterator();
00136                 for (;aNodeIter->more();) {
00137                         const SMDS_MeshNode* aNode = aNodeIter->next();
00138                         current_node.Set(float(aNode->X()),float(aNode->Y()),float(aNode->Z()));
00139                         vertices.push_back(current_node);
00140                 }
00141 
00142                 MeshCore::MeshFacet aFacet;
00143                 aFacet._aulPoints[0] = 0;aFacet._aulPoints[1] = 1;aFacet._aulPoints[2] = 2;
00144                 faces.push_back(aFacet);
00145                 //Fill the Kernel with the temp smesh structure and delete the current containers
00146                 aMesh.Adopt(vertices,faces);
00147                 MeshCore::MeshEigensystem pca(aMesh);
00148                 pca.Evaluate();
00149                 Base::Matrix4D Trafo = pca.Transform();
00150                 /*Let´s transform the input mesh with the PCA Matrix*/
00151                 inputMesh->getFemMeshPtr()->transformGeometry(Trafo);
00152                 //inputMesh->getFemMeshPtr()->getSMesh()->ExportUNV("C:/PCA_alignment.unv");
00153         } PY_CATCH;
00154 
00155         Py_Return;
00156 }
00157 
00158 static PyObject * calcMeshVolume(PyObject *self, PyObject *args)
00159 {
00160         PyObject *input;
00161 
00162         if (!PyArg_ParseTuple(args, "O",&input))
00163                 return NULL;
00164 
00165         PY_TRY {
00166 
00167                 FemMeshPy *inputMesh = static_cast<FemMeshPy*>(input); 
00168 
00169                 SMDS_VolumeIteratorPtr aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator();
00170                 Base::Vector3d a,b,c,a_b_product,temp,temp1;
00171                 double volume =0.0;
00172                 for (;aVolIter->more();) 
00173                 {
00174                         const SMDS_MeshVolume* aVol = aVolIter->next();
00175                         //To make sure that the volume calculation is based on the ABAQUS element convention
00176                         //The following Node mapping from SMESH to ABAQUS is necessary
00177                         //ABAQUS_Node_Number|SMESH_Node_Number
00178                         //0|0
00179                         //1|2
00180                         //2|1
00181                         //3|3
00182                         //4|6
00183                         //5|5
00184                         //6|4
00185                         //7|8
00186                         //8|9
00187                         //9|7
00188                         //The following coordinates of the little pyramids are based on ABAQUS convention and are numbered from
00189                         //1 to 10
00190                         //1,5,8,7
00191                         temp.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z());
00192                         temp1.Set(aVol->GetNode(0)->X(),aVol->GetNode(0)->Y(),aVol->GetNode(0)->Z());
00193                         a = temp - temp1;
00194                         temp.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z());
00195                         temp1.Set(aVol->GetNode(0)->X(),aVol->GetNode(0)->Y(),aVol->GetNode(0)->Z());
00196                         b = temp - temp1;
00197                         temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); 
00198                         temp1.Set(aVol->GetNode(0)->X(),aVol->GetNode(0)->Y(),aVol->GetNode(0)->Z());
00199                         c = temp - temp1;
00200                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00201                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00202                         //5,9,8,7
00203                         temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); 
00204                         temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z());
00205                         a = temp - temp1;
00206                         temp.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z());
00207                         temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z());
00208                         b = temp - temp1;
00209                         temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z());
00210                         temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z());
00211                         c = temp - temp1;
00212                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00213                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00214                         //5,2,9,7
00215                         temp.Set(aVol->GetNode(2)->X(),aVol->GetNode(2)->Y(),aVol->GetNode(2)->Z());
00216                         temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z());
00217                         a = temp - temp1;
00218                         temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); 
00219                         temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z());
00220                         b = temp - temp1;
00221                         temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); 
00222                         temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z());
00223                         c = temp - temp1;
00224                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00225                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00226                         //2,6,9,7
00227                         temp.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z());
00228                         temp1.Set(aVol->GetNode(2)->X(),aVol->GetNode(2)->Y(),aVol->GetNode(2)->Z());
00229                         a = temp - temp1;
00230                         temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); 
00231                         temp1.Set(aVol->GetNode(2)->X(),aVol->GetNode(2)->Y(),aVol->GetNode(2)->Z());
00232                         b = temp - temp1;
00233                         temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z());
00234                         temp1.Set(aVol->GetNode(2)->X(),aVol->GetNode(2)->Y(),aVol->GetNode(2)->Z());
00235                         c = temp - temp1;
00236                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00237                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00238                         //9,6,10,7
00239                         temp.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z()); 
00240                         temp1.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z());
00241                         a = temp - temp1;
00242                         temp.Set(aVol->GetNode(7)->X(),aVol->GetNode(7)->Y(),aVol->GetNode(7)->Z()); 
00243                         temp1.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z());
00244                         b = temp - temp1;
00245                         temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); 
00246                         temp1.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z());
00247                         c = temp - temp1;
00248                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00249                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00250                         //6,3,10,7
00251                         temp.Set(aVol->GetNode(1)->X(),aVol->GetNode(1)->Y(),aVol->GetNode(1)->Z());
00252                         temp1.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z());
00253                         a = temp - temp1;
00254                         temp.Set(aVol->GetNode(7)->X(),aVol->GetNode(7)->Y(),aVol->GetNode(7)->Z()); 
00255                         temp1.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z());
00256                         b = temp - temp1;
00257                         temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); 
00258                         temp1.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z());
00259                         c = temp - temp1;
00260                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00261                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00262                         //8,9,10,7
00263                         temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); 
00264                         temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z());
00265                         a = temp - temp1;
00266                         temp.Set(aVol->GetNode(7)->X(),aVol->GetNode(7)->Y(),aVol->GetNode(7)->Z()); 
00267                         temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z());
00268                         b = temp - temp1;
00269                         temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); 
00270                         temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z());
00271                         c = temp - temp1;
00272                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00273                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00274                         //8,9,10,4
00275                         temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z());
00276                         temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z());
00277                         a = temp - temp1;
00278                         temp.Set(aVol->GetNode(7)->X(),aVol->GetNode(7)->Y(),aVol->GetNode(7)->Z()); 
00279                         temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z());
00280                         b = temp - temp1;
00281                         temp.Set(aVol->GetNode(3)->X(),aVol->GetNode(3)->Y(),aVol->GetNode(3)->Z()); 
00282                         temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z());
00283                         c = temp - temp1;
00284                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00285                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00286                 }
00287                 Py::Float py_volume(volume); 
00288                 return Py::new_reference_to(py_volume);
00289 
00290         } PY_CATCH;
00291 
00292         Py_Return;
00293 }
00294 
00295 static PyObject * checkBB(PyObject *self, PyObject *args)
00296 {
00297         PyObject *input;
00298         PyObject* plm=0;
00299         float billet_thickness;
00300         bool oversize = false;
00301 
00302     if (!PyArg_ParseTuple(args, "O|O!f", &input,&(Base::PlacementPy::Type),&plm,&billet_thickness))
00303         return NULL;
00304 
00305     try {
00306         Base::Placement* placement = 0;
00307         if (plm) {
00308             placement = static_cast<Base::PlacementPy*>(plm)->getPlacementPtr();
00309 
00310         }
00311                 Base::Vector3d current_node;
00312         Base::Matrix4D matrix = placement->toMatrix();
00313         FemMeshPy *inputMesh = static_cast<FemMeshPy*>(input); 
00314                 SMDS_NodeIteratorPtr aNodeIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->nodesIterator();
00315                 for (;aNodeIter->more();) {
00316                         const SMDS_MeshNode* aNode = aNodeIter->next();
00317                         current_node.Set(float(aNode->X()),float(aNode->Y()),float(aNode->Z()));
00318             current_node = matrix * current_node;
00319                         if(current_node.z > billet_thickness || current_node.z < 0.0)
00320                         {
00321                                 //lets jump out of the function as soon as we find a 
00322                                 //Node that is higher or lower than billet thickness
00323                                 oversize = true;
00324                                 Py::Boolean py_oversize(oversize); 
00325                                 return Py::new_reference_to(py_oversize);
00326                         }
00327                 }
00328                 Py::Boolean py_oversize(oversize); 
00329                 return Py::new_reference_to(py_oversize);
00330     }
00331     catch (const std::exception& e) {
00332         PyErr_SetString(PyExc_Exception, e.what());
00333         return 0;
00334     }
00335     Py_Return;
00336 }
00337 
00338 
00339 
00340 
00341 static PyObject * getBoundary_Conditions(PyObject *self, PyObject *args)
00342 {
00343         PyObject *input;
00344         Py::List boundary_nodes;
00345 
00346         if (!PyArg_ParseTuple(args, "O",&input))
00347                 return NULL;
00348 
00349         PY_TRY {
00350 
00351                 FemMeshPy *inputMesh = static_cast<FemMeshPy*>(input); 
00352                 MeshCore::MeshKernel aMesh;
00353                 MeshCore::MeshPointArray vertices;
00354                 vertices.clear();
00355                 MeshCore::MeshFacetArray faces;
00356                 faces.clear();
00357                 MeshCore::MeshPoint current_node;
00358                 SMDS_NodeIteratorPtr aNodeIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->nodesIterator();
00359                 SMDS_VolumeIteratorPtr aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator();
00360                 for (;aNodeIter->more();) {
00361                         const SMDS_MeshNode* aNode = aNodeIter->next();
00362                         current_node.Set(float(aNode->X()),float(aNode->Y()),float(aNode->Z()));
00363                         vertices.push_back(current_node);
00364                 }
00365                 MeshCore::MeshFacet aFacet;
00366                 aFacet._aulPoints[0] = 0;aFacet._aulPoints[1] = 1;aFacet._aulPoints[2] = 2;
00367                 faces.push_back(aFacet);
00368                 //Fill the Kernel with the temp smesh structure and delete the current containers
00369                 aMesh.Adopt(vertices,faces);
00370                 Base::BoundBox3f aBBox;
00371                 aBBox = aMesh.GetBoundBox();
00372 
00373                 float dist_length;
00374         Base::Vector3f dist;
00375                 int minNodeID,maxNodeID,midNodeID;
00376                 dist_length = FLOAT_MAX;
00377                 
00378                 aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator();
00379                 //We only search in non midside nodes (equals the first four nodes of each element)
00380                 for (;aVolIter->more();) 
00381                 {
00382                         const SMDS_MeshVolume* aVol = aVolIter->next();
00383                         for (unsigned j=0;j<4;j++)
00384                         {
00385                                 const SMDS_MeshNode* aNode = aVol->GetNode(j);
00386                                 //Calc distance between the lower left corner and the most next point of the mesh
00387                                 dist.x = float(aNode->X())-aBBox.MinX;dist.y = float(aNode->Y())-aBBox.MinY;dist.z = float(aNode->Z())-aBBox.MinZ;
00388                                 if(dist.Length()<dist_length)
00389                                 {
00390                                         minNodeID = aNode->GetID();
00391                                         dist_length = dist.Length();
00392                                 }
00393                         }
00394                 }
00395 
00396                 boundary_nodes.append(Py::Int(minNodeID));
00397                 
00398                 dist_length = FLOAT_MAX;
00399                 aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator();
00400                 for (;aVolIter->more();) 
00401                 {
00402                         const SMDS_MeshVolume* aVol = aVolIter->next();
00403                         for (unsigned j=0;j<4;j++)
00404                         {
00405                                 const SMDS_MeshNode* aNode = aVol->GetNode(j);
00406                         //Calc distance between the lower right corner and the most next point of the mesh
00407                         dist.x = float(aNode->X())-aBBox.MaxX;dist.y = float(aNode->Y())-aBBox.MinY;dist.z = float(aNode->Z())-aBBox.MinZ;
00408                         if(dist.Length()<dist_length)
00409                         {
00410                                 midNodeID = aNode->GetID();
00411                                 dist_length = dist.Length();
00412                         }
00413                         }
00414                 }
00415                 boundary_nodes.append(Py::Int(midNodeID));
00416                 
00417                 dist_length = FLOAT_MAX;
00418                 aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator();
00419                 for (;aVolIter->more();) 
00420                 {
00421                         const SMDS_MeshVolume* aVol = aVolIter->next();
00422                         for (unsigned j=0;j<4;j++)
00423                         {
00424                                 const SMDS_MeshNode* aNode = aVol->GetNode(j);
00425                         //Calc distance between the lowest lower right corner and the most next point of the mesh
00426                         dist.x = float(aNode->X())-aBBox.MinX;dist.y = float(aNode->Y())-aBBox.MaxY;dist.z = float(aNode->Z())-aBBox.MinZ;
00427                         if(dist.Length()<dist_length)
00428                         {
00429                                 maxNodeID = aNode->GetID();
00430                                 dist_length = dist.Length();
00431                         }
00432                         }
00433                 }
00434                 boundary_nodes.append(Py::Int(maxNodeID));
00435 
00436 
00437                 
00438 
00439                 return Py::new_reference_to(boundary_nodes);
00440                 
00441                 
00442         } PY_CATCH;
00443 
00444         Py_Return;
00445 }
00446 
00447 
00448 
00449 
00450 static PyObject * minBoundingBox(PyObject *self, PyObject *args)
00451 {
00452         
00453         PyObject *input;
00454 
00455         if (!PyArg_ParseTuple(args, "O",&input))
00456                 return NULL;
00457 
00458         PY_TRY {
00459                 FemMeshPy *inputMesh = static_cast<FemMeshPy*>(input); 
00460                 MeshCore::MeshKernel aMesh;
00461                 MeshCore::MeshPointArray vertices;
00462                 vertices.clear();
00463                 MeshCore::MeshFacetArray faces;
00464                 faces.clear();
00465                 MeshCore::MeshPoint current_node;
00466                 SMDS_NodeIteratorPtr aNodeIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->nodesIterator();
00467                 for (;aNodeIter->more();) {
00468                         const SMDS_MeshNode* aNode = aNodeIter->next();
00469                         current_node.Set(float(aNode->X()),float(aNode->Y()),float(aNode->Z()));
00470                         vertices.push_back(current_node);
00471                 }
00472                 MeshCore::MeshFacet aFacet;
00473                 aFacet._aulPoints[0] = 0;aFacet._aulPoints[1] = 1;aFacet._aulPoints[2] = 2;
00474                 faces.push_back(aFacet);
00475                 //Fill the Kernel with the temp smesh structure and delete the current containers
00476                 aMesh.Adopt(vertices,faces);
00477 
00479                 //Now do Monte Carlo to minimize the BBox of the part
00480                 //Use Quaternions for the rotation stuff
00481                 Base::Rotation rotatex,rotatey,rotatez;
00482                 const Base::Vector3d rotate_axis_x(1.0,0.0,0.0),rotate_axis_y(0.0,1.0,0.0),rotate_axis_z(0.0,0.0,1.0);
00483 
00484                 //Rotate around each axes and choose the settings for the min bbox
00485                 Base::Matrix4D final_trafo;
00486                 Base::BoundBox3f aBBox,min_bbox;
00487                 double volumeBBOX,min_volumeBBOX;
00488                 //Get the current min_volumeBBOX and look if we find a lower one
00489                 aBBox = aMesh.GetBoundBox();
00490                 min_volumeBBOX = aBBox.LengthX()*aBBox.LengthY()*aBBox.LengthZ();
00491 
00492                 MeshCore::MeshKernel atempkernel;
00493 
00494                 float it_steps=10.0;
00495                 double step_size;
00496                 double alpha_x=0.0,alpha_y=0.0,alpha_z=0.0;
00497                 double perfect_ax=0.0,perfect_ay=0.0,perfect_az=0.0;
00498 
00499                 //Do a Monte Carlo approach and start from the Principal Axis System
00500                 //and rotate +/- 60° around each axis in a first iteration
00501                 double  angle_range_min_x=-PI/3.0,angle_range_max_x=PI/3.0,
00502                         angle_range_min_y=-PI/3.0,angle_range_max_y=PI/3.0,
00503                         angle_range_min_z=-PI/3.0,angle_range_max_z=PI/3.0;
00504 
00505                 //We rotate until we are 0.1° sure to be in the right position
00506                 for (step_size = (2.0*PI/it_steps);step_size>(2.0*PI/3600.0);step_size=(2.0*PI/it_steps))
00507                 {
00508                         for(alpha_x=angle_range_min_x;alpha_x<angle_range_max_x;alpha_x=alpha_x+step_size)
00509                         {
00510                                 rotatex.setValue(rotate_axis_x,alpha_x);
00511                                 for(alpha_y=angle_range_min_y;alpha_y<angle_range_max_y;alpha_y=alpha_y+step_size)
00512                                 {
00513                                         rotatey.setValue(rotate_axis_y,alpha_y);
00514                                         for(alpha_z=angle_range_min_z;alpha_z<angle_range_max_z;alpha_z=alpha_z+step_size)
00515                                         {
00516                                                 rotatez.setValue(rotate_axis_z,alpha_z);
00517                                                 (rotatex*rotatey*rotatez).getValue(final_trafo);
00518                                                 atempkernel = aMesh;
00519                                                 atempkernel.Transform(final_trafo);
00520                                                 aBBox = atempkernel.GetBoundBox();
00521                                                 volumeBBOX = aBBox.LengthX()*aBBox.LengthY()*aBBox.LengthZ();
00522                                                 if (volumeBBOX < min_volumeBBOX)
00523                                                 { 
00524                                                         min_volumeBBOX=volumeBBOX;
00525                                                         perfect_ax=alpha_x;
00526                                                         perfect_ay=alpha_y;
00527                                                         perfect_az=alpha_z;
00528                                                 }
00529                                         }
00530                                 }
00531                         }
00532                         //We found a better position than the PAS, now lets fine tune this position
00533                         //and search only in the corridor +/- step_size for an even better one
00534                         angle_range_min_x = perfect_ax - step_size;
00535                         angle_range_max_x = perfect_ax + step_size;
00536                         angle_range_min_y = perfect_ay - step_size;
00537                         angle_range_max_y = perfect_ay + step_size;
00538                         angle_range_min_z = perfect_az - step_size;
00539                         angle_range_max_z = perfect_az + step_size;
00540                         it_steps = it_steps*float(5.0);
00541                 }
00542 
00544                 //Free Memory
00545                 atempkernel.Clear();
00546                 //Transform the mesh to the evaluated perfect position right now
00547                 rotatex.setValue(rotate_axis_x,perfect_ax);
00548                 rotatey.setValue(rotate_axis_y,perfect_ay);
00549                 rotatez.setValue(rotate_axis_z,perfect_az);
00550                 (rotatex*rotatey*rotatez).getValue(final_trafo);
00551                 aMesh.Transform(final_trafo);
00552                 inputMesh->getFemMeshPtr()->transformGeometry(final_trafo);
00554                 //Now lets also do the movement to the 1st Quadrant in this function
00555                 aBBox = aMesh.GetBoundBox();
00556                 //Get Distance vector from current lower left corner of BBox to origin
00557                 Base::Vector3f dist_vector;
00558                 dist_vector.x = -aBBox.MinX;dist_vector.y=-aBBox.MinY;dist_vector.z=-aBBox.MinZ;
00559                 Base::Matrix4D trans_matrix(
00560                         float(1.0),float(0.0),float(0.0),dist_vector.x,
00561                         float(0.0),float(1.0),float(0.0),dist_vector.y,
00562                         float(0.0),float(0.0),float(1.0),dist_vector.z,
00563                         float(0.0),float(0.0),float(0.0),float(1.0));
00564                 inputMesh->getFemMeshPtr()->transformGeometry(trans_matrix);
00565                 
00566                 //inputMesh->getFemMeshPtr()->getSMesh()->ExportUNV("C:/fine_tuning.unv");
00567 
00568         } PY_CATCH;
00569 
00570         Py_Return;
00571 }
00572 
00573 
00574 static PyObject * importer(PyObject *self, PyObject *args)
00575 {
00576     const char* Name;
00577     const char* DocName=0;
00578 
00579     if (!PyArg_ParseTuple(args, "s|s",&Name,&DocName))
00580         return NULL;
00581 
00582     PY_TRY {
00583         App::Document *pcDoc = 0;
00584         if (DocName)
00585             pcDoc = App::GetApplication().getDocument(DocName);
00586         else
00587             pcDoc = App::GetApplication().getActiveDocument();
00588 
00589         if (!pcDoc) {
00590             pcDoc = App::GetApplication().newDocument(DocName);
00591         }
00592 
00593         std::auto_ptr<FemMesh> mesh(new FemMesh);
00594         mesh->read(Name);
00595         Base::FileInfo file(Name);
00596         FemMeshObject *pcFeature = static_cast<FemMeshObject *>
00597             (pcDoc->addObject("Fem::FemMeshObject", file.fileNamePure().c_str()));
00598         pcFeature->Label.setValue(file.fileNamePure().c_str());
00599         pcFeature->FemMesh.setValuePtr(mesh.get());
00600         (void)mesh.release();
00601         pcFeature->purgeTouched();
00602     } PY_CATCH;
00603 
00604     Py_Return;
00605 }
00606 
00607 
00608 static PyObject * import_NASTRAN(PyObject *self, PyObject *args)
00609 {
00610         const char* filename_input, *filename_output;
00611         if (!PyArg_ParseTuple(args, "ss",&filename_input,&filename_output))
00612                 return NULL;
00613 
00614         PY_TRY {
00615 
00616                 std::ifstream inputfile;
00617 
00618                 //Für Debugoutput
00619                 ofstream anOutput;
00620                 anOutput.open("c:/time_measurement.txt");
00621                 time_t seconds1,seconds2;
00622                 
00623                 inputfile.open(filename_input);
00624                 if (!inputfile.is_open())  //Exists...?
00625                 {
00626                         std::cerr << "File not found. Exiting..." << std::endl;
00627                         return NULL;
00628                 }
00629 
00630                 //Return the line pointer to the beginning of the file
00631                 inputfile.seekg(std::ifstream::beg);
00632                 std::string line1,line2,temp;
00633                 std::stringstream astream;
00634                 std::vector<unsigned int> tetra_element;
00635                 std::vector<unsigned int> element_id;
00636                 element_id.clear();
00637                 std::vector<std::vector<unsigned int> > all_elements;
00638                 std::vector<std::vector<unsigned int> >::iterator all_element_iterator;
00639                 std::vector<unsigned int>::iterator one_element_iterator;
00640                 all_elements.clear();
00641                 MeshCore::MeshKernel aMesh;
00642                 MeshCore::MeshPointArray vertices;
00643                 vertices.clear();
00644                 MeshCore::MeshFacetArray faces;
00645                 faces.clear();
00646                 MeshCore::MeshPoint current_node;
00647 
00648                 seconds1 = time(NULL);
00649                 do
00650                 {
00651                         std::getline(inputfile,line1);
00652                         if (line1.size() == 0) continue;
00653                         if (line1.find("GRID*")!= std::string::npos) //We found a Grid line
00654                         {
00655                                 //Now lets extract the GRID Points = Nodes
00656                                 //As each GRID Line consists of two subsequent lines we have to
00657                                 //take care of that as well
00658                                 std::getline(inputfile,line2);
00659                                 //Extract X Value
00660                                 current_node.x = float(atof(line1.substr(40,56).c_str()));
00661                                 //Extract Y Value
00662                                 current_node.y = float(atof(line1.substr(56,72).c_str()));
00663                                 //Extract Z Value
00664                                 current_node.z = float(atof(line2.substr(8,24).c_str()));
00665                                 
00666                                 vertices.push_back(current_node);
00667                         }
00668                         else if (line1.find("CTETRA")!= std::string::npos)
00669                         {
00670                                 tetra_element.clear();
00671                                 //Lets extract the elements
00672                                 //As each Element Line consists of two subsequent lines as well 
00673                                 //we have to take care of that
00674                                 //At a first step we only extract Quadratic Tetrahedral Elements
00675                                 std::getline(inputfile,line2);
00676                                 element_id.push_back(atoi(line1.substr(8,16).c_str()));
00677                                 tetra_element.push_back(atoi(line1.substr(24,32).c_str()));
00678                                 tetra_element.push_back(atoi(line1.substr(32,40).c_str()));
00679                                 tetra_element.push_back(atoi(line1.substr(40,48).c_str()));
00680                                 tetra_element.push_back(atoi(line1.substr(48,56).c_str()));
00681                                 tetra_element.push_back(atoi(line1.substr(56,64).c_str()));
00682                                 tetra_element.push_back(atoi(line1.substr(64,72).c_str()));
00683                                 tetra_element.push_back(atoi(line2.substr(8,16).c_str()));
00684                                 tetra_element.push_back(atoi(line2.substr(16,24).c_str()));
00685                                 tetra_element.push_back(atoi(line2.substr(24,32).c_str()));
00686                                 tetra_element.push_back(atoi(line2.substr(32,40).c_str()));
00687                         
00688                                 all_elements.push_back(tetra_element);
00689                         }
00690                         
00691                 }
00692                 while (inputfile.good());
00693                 inputfile.close();
00694 
00695                 seconds2 = time(NULL);
00696 
00697                 anOutput << seconds2-seconds1 <<" for Parsing the input file"<<std::endl;
00698                 //Now lets perform some minimization routines to min the bbox volume
00699                 //To evaluate the COG and principal axes we have to generate a "Mesh" out of the points
00700                 //therefore at least one face is required
00701                 MeshCore::MeshFacet aFacet;
00702                 aFacet._aulPoints[0] = 0;aFacet._aulPoints[1] = 1;aFacet._aulPoints[2] = 2;
00703                 faces.push_back(aFacet);
00704                 //Fill the Kernel with the structure
00705             aMesh.Assign(vertices,faces);
00706                 
00707                 seconds1 = time(NULL);
00708                 //First lets get the COG and the principal axes and transfer the mesh there
00709                 MeshCore::MeshEigensystem pca(aMesh);
00710                 pca.Evaluate();
00711                 Base::Matrix4D Trafo =  pca.Transform();
00712                 aMesh.Transform(Trafo);
00714                 //To get the center of gravity we have to build the inverse of the pca Matrix
00715                 pca.Evaluate();
00716                 Trafo =  pca.Transform();
00717                 Trafo.inverse();
00718                 Base::Vector3f cog;
00719                 cog.x = float(Trafo[0][3]);cog.y = float(Trafo[1][3]);cog.z = float(Trafo[2][3]);
00721                 //Now do Monte Carlo to minimize the BBox of the part
00722                 //Use Quaternions for the rotation stuff
00723                 Base::Rotation rotatex,rotatey,rotatez;
00724                 const Base::Vector3d rotate_axis_x(1.0,0.0,0.0),rotate_axis_y(0.0,1.0,0.0),rotate_axis_z(0.0,0.0,1.0);
00725                 
00726                 //Rotate around each axes and choose the settings for the min bbox
00727                 Base::Matrix4D final_trafo;
00728                 Base::BoundBox3f aBBox,min_bbox;
00729 
00730 
00731                 double volumeBBOX,min_volumeBBOX;
00732 
00733                 //Get the current min_volumeBBOX and look if we find a lower one
00734                 aBBox = aMesh.GetBoundBox();
00735                 min_volumeBBOX = aBBox.LengthX()*aBBox.LengthY()*aBBox.LengthZ();
00736 
00737         
00738                 MeshCore::MeshKernel atempkernel;
00739 
00740                 float it_steps=10.0;
00741                 double step_size;
00742                 double alpha_x=0.0,alpha_y=0.0,alpha_z=0.0;
00743                 double perfect_ax=0.0,perfect_ay=0.0,perfect_az=0.0;
00744 
00745                 //Do a Monte Carlo approach and start from the Principal Axis System
00746                 //and rotate +/- 60° around each axis in a first iteration
00747                 double  angle_range_min_x=-PI/3.0,angle_range_max_x=PI/3.0,
00748                                 angle_range_min_y=-PI/3.0,angle_range_max_y=PI/3.0,
00749                                 angle_range_min_z=-PI/3.0,angle_range_max_z=PI/3.0;
00750                 
00751                 for (step_size = (2.0*PI/it_steps);step_size>(2.0*PI/360.0);step_size=(2.0*PI/it_steps))
00752                 {
00753                         for(alpha_x=angle_range_min_x;alpha_x<angle_range_max_x;alpha_x=alpha_x+step_size)
00754                         {
00755                                 rotatex.setValue(rotate_axis_x,alpha_x);
00756                                 for(alpha_y=angle_range_min_y;alpha_y<angle_range_max_y;alpha_y=alpha_y+step_size)
00757                                 {
00758                                         rotatey.setValue(rotate_axis_y,alpha_y);
00759                                         for(alpha_z=angle_range_min_z;alpha_z<angle_range_max_z;alpha_z=alpha_z+step_size)
00760                                         {
00761                                                 rotatez.setValue(rotate_axis_z,alpha_z);
00762                                                 (rotatex*rotatey*rotatez).getValue(final_trafo);
00763                                                 atempkernel = aMesh;
00764                                                 atempkernel.Transform(final_trafo);
00765                                                 aBBox = atempkernel.GetBoundBox();
00766                                                 volumeBBOX = aBBox.LengthX()*aBBox.LengthY()*aBBox.LengthZ();
00767                                                 if (volumeBBOX < min_volumeBBOX)
00768                                                 { 
00769                                                         min_volumeBBOX=volumeBBOX;
00770                                                         perfect_ax=alpha_x;
00771                                                         perfect_ay=alpha_y;
00772                                                         perfect_az=alpha_z;
00773                                                 }
00774                                         }
00775                                 }
00776                         }
00777                         //We found a better position than the PAS, now lets fine tune this position
00778                         //and search only in the corridor +/- step_size for an even better one
00779                         angle_range_min_x = perfect_ax - step_size;
00780                         angle_range_max_x = perfect_ax + step_size;
00781                         angle_range_min_y = perfect_ay - step_size;
00782                         angle_range_max_y = perfect_ay + step_size;
00783                         angle_range_min_z = perfect_az - step_size;
00784                         angle_range_max_z = perfect_az + step_size;
00785                         it_steps = it_steps*float(5.0);
00786                 }
00787 
00789                 
00790                 //Free Memory
00791                 atempkernel.Clear();
00792 
00793                 //Transform the mesh to the evaluated perfect position right now
00794                 rotatex.setValue(rotate_axis_x,perfect_ax);
00795                 rotatey.setValue(rotate_axis_y,perfect_ay);
00796                 rotatez.setValue(rotate_axis_z,perfect_az);
00797                 (rotatex*rotatey*rotatez).getValue(final_trafo);
00798                 aMesh.Transform(final_trafo);
00799 
00800                 anOutput << "perfect angles " << perfect_ax << "," << perfect_ay << "," << perfect_az << std::endl;
00801 
00802 
00803                 //Move Mesh to stay fully in the positive octant
00804                 aBBox = aMesh.GetBoundBox();
00805                 //Get Distance vector from current lower left corner of BBox to origin
00806                 Base::Vector3f dist_vector;
00807                 dist_vector.x = -aBBox.MinX;dist_vector.y=-aBBox.MinY;dist_vector.z=-aBBox.MinZ;
00808                 Base::Matrix4D trans_matrix(
00809                         float(1.0),float(0.0),float(0.0),dist_vector.x,
00810                         float(0.0),float(1.0),float(0.0),dist_vector.y,
00811                         float(0.0),float(0.0),float(1.0),dist_vector.z,
00812                         float(0.0),float(0.0),float(0.0),float(1.0));
00813                 aMesh.Transform(trans_matrix);
00815                 seconds2=time(NULL);
00816                 anOutput << seconds2-seconds1 << " seconds for the min bounding box stuff" << std::endl;
00817                 
00818                 seconds1 = time(NULL);
00819                 //Try to build up an own mesh kernel within SMESH to export the 
00820                 //perfect positioned mesh right now
00821                 SMESH_Gen* meshgen = new SMESH_Gen();
00822                 SMESH_Mesh* mesh = meshgen->CreateMesh(1, true);
00823                 SMESHDS_Mesh* meshds = mesh->GetMeshDS();
00824                 
00825                 MeshCore::MeshPointIterator anodeiterator(aMesh);
00826                 int j=1;
00827                 for(anodeiterator.Begin(); anodeiterator.More(); anodeiterator.Next())
00828                 {
00829                         meshds->AddNodeWithID((*anodeiterator).x,(*anodeiterator).y,(*anodeiterator).z,j);
00830                         j++;
00831                 }
00832 
00833                 for(unsigned int i=0;i<all_elements.size();i++)
00834                 {
00835                         //Die Reihenfolge wie hier die Elemente hinzugefügt werden ist sehr wichtig. 
00836                         //Ansonsten ist eine konsistente Datenstruktur nicht möglich
00837                         meshds->AddVolumeWithID(
00838                                 meshds->FindNode(all_elements[i][0]),
00839                                 meshds->FindNode(all_elements[i][2]),
00840                                 meshds->FindNode(all_elements[i][1]),
00841                                 meshds->FindNode(all_elements[i][3]),
00842                                 meshds->FindNode(all_elements[i][6]),
00843                                 meshds->FindNode(all_elements[i][5]),
00844                                 meshds->FindNode(all_elements[i][4]),
00845                                 meshds->FindNode(all_elements[i][9]),
00846                                 meshds->FindNode(all_elements[i][7]),
00847                                 meshds->FindNode(all_elements[i][8]),
00848                                 element_id[i]
00849                         );
00850                 }
00851 
00852                 mesh->ExportUNV(filename_output);
00854                 seconds2 = time(NULL);
00855                 anOutput << seconds2-seconds1 << " seconds for the Mesh Export" << std::endl;
00856 
00857                 //Output also to ABAQUS Input Format
00858                 ofstream anABAQUS_Output;
00859                 anABAQUS_Output.open("d:/abaqus_output.inp");
00860                 anABAQUS_Output << "*Node , NSET=Nall" << std::endl;
00861                 j=1;
00862                 for(anodeiterator.Begin(); anodeiterator.More(); anodeiterator.Next())
00863                 {
00864                         anABAQUS_Output << j <<","
00865                    <<(*anodeiterator).x << "," 
00866                    <<(*anodeiterator).y << ","
00867                    <<(*anodeiterator).z << std::endl;
00868                         j++;
00869                 }
00870                 anABAQUS_Output << "*Element, TYPE=C3D10, ELSET=Eall" << std::endl;
00871                 j=1;
00872                 for(unsigned int i=0;i<all_elements.size();i++)
00873                 {
00874                         //In ABAQUS input format a maximum of 15 Nodes per line is allowed
00875                         anABAQUS_Output 
00876                         <<j <<","
00877                         <<all_elements[i][0]<<","
00878                         <<all_elements[i][1]<<","
00879                         <<all_elements[i][2]<<","
00880                         <<all_elements[i][3]<<","
00881                         <<all_elements[i][4]<<","
00882                         <<all_elements[i][5]<<","
00883                         <<all_elements[i][6]<<","
00884                         <<all_elements[i][7]<<","
00885                         <<all_elements[i][8]<<","
00886                         <<all_elements[i][9]<<std::endl;
00887                         j++;
00888                 }
00889                 anABAQUS_Output.close();
00891 
00892 
00893                 //Calculate Mesh Volume
00894                 //For an accurate Volume Calculation of a quadratic Tetrahedron
00895                 //we have to calculate the Volume of 8 Sub-Tetrahedrons
00896                 Base::Vector3f a,b,c,a_b_product;
00897                 double volume = 0.0;
00898                 seconds1 = time(NULL);
00899                 //Calc Volume with 1/6 * |(a x b) * c|
00900                 for(all_element_iterator = all_elements.begin();all_element_iterator != all_elements.end();all_element_iterator++)
00901                 {
00902                         //1,5,8,7
00903                         a = vertices[all_element_iterator->at(4)-1]-vertices[all_element_iterator->at(0)-1];
00904                         b = vertices[all_element_iterator->at(7)-1]-vertices[all_element_iterator->at(0)-1];
00905                         c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(0)-1];
00906                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00907                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00908                         //5,9,8,7
00909                         a = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(4)-1];
00910                         b = vertices[all_element_iterator->at(7)-1]-vertices[all_element_iterator->at(4)-1];
00911                         c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(4)-1];
00912                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00913                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00914                         //5,2,9,7
00915                         a = vertices[all_element_iterator->at(1)-1]-vertices[all_element_iterator->at(4)-1];
00916                         b = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(4)-1];
00917                         c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(4)-1];
00918                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00919                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00920                         //2,6,9,7
00921                         a = vertices[all_element_iterator->at(5)-1]-vertices[all_element_iterator->at(1)-1];
00922                         b = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(1)-1];
00923                         c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(1)-1];
00924                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00925                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00926                         //9,6,10,7
00927                         a = vertices[all_element_iterator->at(5)-1]-vertices[all_element_iterator->at(8)-1];
00928                         b = vertices[all_element_iterator->at(9)-1]-vertices[all_element_iterator->at(8)-1];
00929                         c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(8)-1];
00930                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00931                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00932                         //6,3,10,7
00933                         a = vertices[all_element_iterator->at(2)-1]-vertices[all_element_iterator->at(5)-1];
00934                         b = vertices[all_element_iterator->at(9)-1]-vertices[all_element_iterator->at(5)-1];
00935                         c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(5)-1];
00936                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00937                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00938                         //8,9,10,7
00939                         a = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(7)-1];
00940                         b = vertices[all_element_iterator->at(9)-1]-vertices[all_element_iterator->at(7)-1];
00941                         c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(7)-1];
00942                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00943                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00944                         //8,9,10,4
00945                         a = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(7)-1];
00946                         b = vertices[all_element_iterator->at(9)-1]-vertices[all_element_iterator->at(7)-1];
00947                         c = vertices[all_element_iterator->at(3)-1]-vertices[all_element_iterator->at(7)-1];
00948                         a_b_product.x = a.y*b.z-b.y*a.z;a_b_product.y = a.z*b.x-b.z*a.x;a_b_product.z = a.x*b.y-b.x*a.y;
00949                         volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z));
00950                 
00951                 }
00952                 seconds2=time(NULL);
00953                 anOutput << seconds2-seconds1 << " seconds for Volume Calculation  " << "Volumen " << volume/1000.0/1000.0/1000.0 << " in m^3" << std::endl;
00954                 anOutput  << "Volumen der BBox" << min_volumeBBOX/1000.0/1000.0/1000.0 << std::endl;
00955                 anOutput << "Fly to Buy Ratio: " << min_volumeBBOX / volume << std::endl;
00956                 anOutput.close();
00957                 
00958         } PY_CATCH;
00959 
00960         Py_Return;
00961 }
00962 
00963 
00964 static PyObject * exporter(PyObject *self, PyObject *args)
00965 {
00966     PyObject* object;
00967     const char* filename;
00968     if (!PyArg_ParseTuple(args, "Os",&object,&filename))
00969         return NULL;
00970 
00971     PY_TRY {
00972         Py::List list(object);
00973         Base::Type meshId = Base::Type::fromName("Fem::FemMeshObject");
00974         for (Py::List::iterator it = list.begin(); it != list.end(); ++it) {
00975             PyObject* item = (*it).ptr();
00976             if (PyObject_TypeCheck(item, &(App::DocumentObjectPy::Type))) {
00977                 App::DocumentObject* obj = static_cast<App::DocumentObjectPy*>(item)->getDocumentObjectPtr();
00978                 if (obj->getTypeId().isDerivedFrom(meshId)) {
00979                     static_cast<FemMeshObject*>(obj)->FemMesh.getValue().write(filename);
00980                     break;
00981                 }
00982             }
00983         }
00984     } PY_CATCH;
00985 
00986     Py_Return;
00987 }
00988 
00989 // ----------------------------------------------------------------------------
00990 
00991 PyDoc_STRVAR(open_doc,
00992 "open(string) -- Create a new document and a Mesh::Import feature to load the file into the document.");
00993 
00994 PyDoc_STRVAR(inst_doc,
00995 "insert(string|mesh,[string]) -- Load or insert a mesh into the given or active document.");
00996 
00997 PyDoc_STRVAR(export_doc,
00998 "export(list,string) -- Export a list of objects into a single file.");
00999 
01000 /* registration table  */
01001 struct PyMethodDef Fem_methods[] = {
01002     {"open"       ,open ,       METH_VARARGS, open_doc},
01003     {"insert"     ,importer,    METH_VARARGS, inst_doc},
01004     {"export"     ,exporter,    METH_VARARGS, export_doc},
01005     {"read"       ,read,        Py_NEWARGS,   "Read a mesh from a file and returns a Mesh object."},
01006 {"calcMeshVolume", calcMeshVolume, Py_NEWARGS, "Calculate Mesh Volume for C3D10"},      
01007 {"getBoundary_Conditions" , getBoundary_Conditions, Py_NEWARGS, "Get Boundary Conditions for Residual Stress Calculation"},
01008         {"SMESH_PCA" , SMESH_PCA, Py_NEWARGS, "Get a Matrix4D related to the PCA of a Mesh Object"},
01009         {"import_NASTRAN",import_NASTRAN, Py_NEWARGS, "Test"},
01010         {"minBoundingBox",minBoundingBox,Py_NEWARGS,"Minimize the Bounding Box and reorient the mesh to the 1st Quadrant"},
01011         {"checkBB",checkBB,Py_NEWARGS,"Check if the nodal z-values are still in the prescribed range"},
01012     {NULL, NULL}  /* sentinel */
01013 };
Generated on Wed Nov 23 18:59:55 2011 for FreeCAD by  doxygen 1.6.1