diff --git a/src/Mod/Fem/App/AppFemPy.cpp b/src/Mod/Fem/App/AppFemPy.cpp index b29d57097..36cd37e87 100755 --- a/src/Mod/Fem/App/AppFemPy.cpp +++ b/src/Mod/Fem/App/AppFemPy.cpp @@ -35,9 +35,9 @@ #include #include #include -#include -#include -#include +//#include +//#include +//#include #include #include @@ -140,535 +140,535 @@ show(PyObject *self, PyObject *args) } -static PyObject * SMESH_PCA(PyObject *self, PyObject *args) -{ - PyObject *input; - - if (!PyArg_ParseTuple(args, "O",&input)) - return NULL; - - PY_TRY { - - FemMeshPy *inputMesh = static_cast(input); - MeshCore::MeshKernel aMesh; - MeshCore::MeshPointArray vertices; - vertices.clear(); - MeshCore::MeshFacetArray faces; - faces.clear(); - MeshCore::MeshPoint current_node; - SMDS_NodeIteratorPtr aNodeIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->nodesIterator(); - for (;aNodeIter->more();) { - const SMDS_MeshNode* aNode = aNodeIter->next(); - current_node.Set(float(aNode->X()),float(aNode->Y()),float(aNode->Z())); - vertices.push_back(current_node); - } - - MeshCore::MeshFacet aFacet; - aFacet._aulPoints[0] = 0;aFacet._aulPoints[1] = 1;aFacet._aulPoints[2] = 2; - faces.push_back(aFacet); - //Fill the Kernel with the temp smesh structure and delete the current containers - aMesh.Adopt(vertices,faces); - MeshCore::MeshEigensystem pca(aMesh); - pca.Evaluate(); - Base::Matrix4D Trafo = pca.Transform(); - /*Let´s transform the input mesh with the PCA Matrix*/ - inputMesh->getFemMeshPtr()->transformGeometry(Trafo); - //inputMesh->getFemMeshPtr()->getSMesh()->ExportUNV("C:/Temp/PCA_alignment.unv"); - //Now lets check if the smallest dimension of the BBox is oriented towards the Z-Axis. If not, lets rotate it around the X or Y axis - //Use the SMESH structure for that - // aMesh.Transform(Trafo); - - //Base::Rotation rotatex,rotatey,rotatez; - //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); - //double bbox_length_x,bbox_length_y,bbox_length_z; - ////Rotate around the each axes and choose the settings for the min bbox - //Base::Matrix4D final_trafo; - //Base::BoundBox3f aBBox; - ////Get the current BBOX and look for the size - //aBBox = aMesh.GetBoundBox(); - //bbox_length_x = aBBox.LengthX();bbox_length_y = aBBox.LengthY();bbox_length_z = aBBox.LengthZ(); - ////Now do the rotation stuff - //if (bbox_length_z < bbox_length_x && bbox_length_z < bbox_length_y) - // Py_Return; - //else if ( - - - //MeshCore::MeshKernel atempkernel; - - //float it_steps=10.0; - //double step_size; - //double alpha_x=0.0,alpha_y=0.0,alpha_z=0.0; - //double perfect_ax=0.0,perfect_ay=0.0,perfect_az=0.0; - - ////Do a Monte Carlo approach and start from the Principal Axis System - ////and rotate +/- 60° around each axis in a first iteration - //double angle_range_min_x=-PI/3.0,angle_range_max_x=PI/3.0, - // angle_range_min_y=-PI/3.0,angle_range_max_y=PI/3.0, - // angle_range_min_z=-PI/3.0,angle_range_max_z=PI/3.0; - - ////We rotate until we are 0.1° sure to be in the right position - //for (step_size = (2.0*PI/it_steps);step_size>(2.0*PI/3600.0);step_size=(2.0*PI/it_steps)) - //{ - // for(alpha_x=angle_range_min_x;alpha_x(input); - - SMDS_VolumeIteratorPtr aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator(); - Base::Vector3d a,b,c,a_b_product,temp,temp1; - double volume =0.0; - for (;aVolIter->more();) - { - const SMDS_MeshVolume* aVol = aVolIter->next(); - //To make sure that the volume calculation is based on the ABAQUS element convention - //The following Node mapping from SMESH to ABAQUS is necessary - //ABAQUS_Node_Number|SMESH_Node_Number - //0|0 - //1|2 - //2|1 - //3|3 - //4|6 - //5|5 - //6|4 - //7|8 - //8|9 - //9|7 - //The following coordinates of the little pyramids are based on ABAQUS convention and are numbered from - //1 to 10 - //1,5,8,7 - temp.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); - temp1.Set(aVol->GetNode(0)->X(),aVol->GetNode(0)->Y(),aVol->GetNode(0)->Z()); - a = temp - temp1; - temp.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); - temp1.Set(aVol->GetNode(0)->X(),aVol->GetNode(0)->Y(),aVol->GetNode(0)->Z()); - b = temp - temp1; - temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); - temp1.Set(aVol->GetNode(0)->X(),aVol->GetNode(0)->Y(),aVol->GetNode(0)->Z()); - c = temp - temp1; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - //5,9,8,7 - temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); - temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); - a = temp - temp1; - temp.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); - temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); - b = temp - temp1; - temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); - temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); - c = temp - temp1; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - //5,2,9,7 - temp.Set(aVol->GetNode(2)->X(),aVol->GetNode(2)->Y(),aVol->GetNode(2)->Z()); - temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); - a = temp - temp1; - temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); - temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); - b = temp - temp1; - temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); - temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); - c = temp - temp1; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - //2,6,9,7 - temp.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z()); - temp1.Set(aVol->GetNode(2)->X(),aVol->GetNode(2)->Y(),aVol->GetNode(2)->Z()); - a = temp - temp1; - temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); - temp1.Set(aVol->GetNode(2)->X(),aVol->GetNode(2)->Y(),aVol->GetNode(2)->Z()); - b = temp - temp1; - temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); - temp1.Set(aVol->GetNode(2)->X(),aVol->GetNode(2)->Y(),aVol->GetNode(2)->Z()); - c = temp - temp1; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - //9,6,10,7 - temp.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z()); - temp1.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); - a = temp - temp1; - temp.Set(aVol->GetNode(7)->X(),aVol->GetNode(7)->Y(),aVol->GetNode(7)->Z()); - temp1.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); - b = temp - temp1; - temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); - temp1.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); - c = temp - temp1; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - //6,3,10,7 - temp.Set(aVol->GetNode(1)->X(),aVol->GetNode(1)->Y(),aVol->GetNode(1)->Z()); - temp1.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z()); - a = temp - temp1; - temp.Set(aVol->GetNode(7)->X(),aVol->GetNode(7)->Y(),aVol->GetNode(7)->Z()); - temp1.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z()); - b = temp - temp1; - temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); - temp1.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z()); - c = temp - temp1; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - //8,9,10,7 - temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); - temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); - a = temp - temp1; - temp.Set(aVol->GetNode(7)->X(),aVol->GetNode(7)->Y(),aVol->GetNode(7)->Z()); - temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); - b = temp - temp1; - temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); - temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); - c = temp - temp1; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - //8,9,10,4 - temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); - temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); - a = temp - temp1; - temp.Set(aVol->GetNode(7)->X(),aVol->GetNode(7)->Y(),aVol->GetNode(7)->Z()); - temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); - b = temp - temp1; - temp.Set(aVol->GetNode(3)->X(),aVol->GetNode(3)->Y(),aVol->GetNode(3)->Z()); - temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); - c = temp - temp1; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - } - Py::Float py_volume(volume); - return Py::new_reference_to(py_volume); - - } PY_CATCH; - - Py_Return; -} - -static PyObject * checkBB(PyObject *self, PyObject *args) -{ - PyObject *input; - PyObject* plm=0; - float billet_thickness; - bool oversize = false; - - if (!PyArg_ParseTuple(args, "O|O!f", &input,&(Base::PlacementPy::Type),&plm,&billet_thickness)) - return NULL; - - try { - Base::Placement* placement = 0; - if (plm) { - placement = static_cast(plm)->getPlacementPtr(); - - } - Base::Vector3d current_node; - Base::Matrix4D matrix = placement->toMatrix(); - FemMeshPy *inputMesh = static_cast(input); - SMDS_NodeIteratorPtr aNodeIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->nodesIterator(); - for (;aNodeIter->more();) { - const SMDS_MeshNode* aNode = aNodeIter->next(); - current_node.Set(float(aNode->X()),float(aNode->Y()),float(aNode->Z())); - current_node = matrix * current_node; - if(current_node.z > billet_thickness || current_node.z < -0.1) - { - //lets jump out of the function as soon as we find a - //Node that is higher or lower than billet thickness - oversize = true; - Py::Boolean py_oversize(oversize); - return Py::new_reference_to(py_oversize); - } - } - Py::Boolean py_oversize(oversize); - return Py::new_reference_to(py_oversize); - } - catch (const std::exception& e) { - PyErr_SetString(PyExc_Exception, e.what()); - return 0; - } - Py_Return; -} - - - - -static PyObject * getBoundary_Conditions(PyObject *self, PyObject *args) -{ - PyObject *input; - Py::List boundary_nodes; - - if (!PyArg_ParseTuple(args, "O!", &(FemMeshPy::Type), &input)) - //if (!PyArg_ParseTuple(args, "O",&input)) - return NULL; - - PY_TRY { - - FemMeshPy *inputMesh = static_cast(input); - MeshCore::MeshKernel aMesh; - MeshCore::MeshPointArray vertices; - vertices.clear(); - MeshCore::MeshFacetArray faces; - faces.clear(); - MeshCore::MeshPoint current_node; - SMDS_NodeIteratorPtr aNodeIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->nodesIterator(); - SMDS_VolumeIteratorPtr aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator(); - for (;aNodeIter->more();) { - const SMDS_MeshNode* aNode = aNodeIter->next(); - current_node.Set(float(aNode->X()),float(aNode->Y()),float(aNode->Z())); - vertices.push_back(current_node); - } - MeshCore::MeshFacet aFacet; - aFacet._aulPoints[0] = 0;aFacet._aulPoints[1] = 1;aFacet._aulPoints[2] = 2; - faces.push_back(aFacet); - //Fill the Kernel with the temp smesh structure and delete the current containers - aMesh.Adopt(vertices,faces); - Base::BoundBox3f aBBox; - aBBox = aMesh.GetBoundBox(); - - float dist_length; - Base::Vector3f dist; - int minNodeID,maxNodeID,midNodeID; - dist_length = FLOAT_MAX; - - aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator(); - //We only search in non midside nodes (equals the first four nodes of each element) - for (;aVolIter->more();) - { - const SMDS_MeshVolume* aVol = aVolIter->next(); - for (unsigned j=0;j<4;j++) - { - const SMDS_MeshNode* aNode = aVol->GetNode(j); - //Calc distance between the lower left corner and the most next point of the mesh - dist.x = float(aNode->X())-aBBox.MinX;dist.y = float(aNode->Y())-aBBox.MinY;dist.z = float(aNode->Z())-aBBox.MinZ; - if(dist.Length()GetID(); - dist_length = dist.Length(); - } - } - } - - boundary_nodes.append(Py::Int(minNodeID)); - - dist_length = FLOAT_MAX; - aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator(); - for (;aVolIter->more();) - { - const SMDS_MeshVolume* aVol = aVolIter->next(); - for (unsigned j=0;j<4;j++) - { - const SMDS_MeshNode* aNode = aVol->GetNode(j); - //Calc distance between the lower right corner and the most next point of the mesh - dist.x = float(aNode->X())-aBBox.MaxX;dist.y = float(aNode->Y())-aBBox.MinY;dist.z = float(aNode->Z())-aBBox.MinZ; - if(dist.Length()GetID(); - dist_length = dist.Length(); - } - } - } - boundary_nodes.append(Py::Int(midNodeID)); - - dist_length = FLOAT_MAX; - aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator(); - for (;aVolIter->more();) - { - const SMDS_MeshVolume* aVol = aVolIter->next(); - for (unsigned j=0;j<4;j++) - { - const SMDS_MeshNode* aNode = aVol->GetNode(j); - //Calc distance between the lowest lower right corner and the most next point of the mesh - dist.x = float(aNode->X())-aBBox.MinX;dist.y = float(aNode->Y())-aBBox.MaxY;dist.z = float(aNode->Z())-aBBox.MinZ; - if(dist.Length()GetID(); - dist_length = dist.Length(); - } - } - } - boundary_nodes.append(Py::Int(maxNodeID)); - - - - - return Py::new_reference_to(boundary_nodes); - - - } PY_CATCH; - - Py_Return; -} - - - - -static PyObject * minBoundingBox(PyObject *self, PyObject *args) -{ - - PyObject *input; - - if (!PyArg_ParseTuple(args, "O",&input)) - return NULL; - - PY_TRY { - FemMeshPy *inputMesh = static_cast(input); - MeshCore::MeshKernel aMesh; - MeshCore::MeshPointArray vertices; - vertices.clear(); - MeshCore::MeshFacetArray faces; - faces.clear(); - MeshCore::MeshPoint current_node; - SMDS_NodeIteratorPtr aNodeIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->nodesIterator(); - for (;aNodeIter->more();) { - const SMDS_MeshNode* aNode = aNodeIter->next(); - current_node.Set(float(aNode->X()),float(aNode->Y()),float(aNode->Z())); - vertices.push_back(current_node); - } - MeshCore::MeshFacet aFacet; - aFacet._aulPoints[0] = 0;aFacet._aulPoints[1] = 1;aFacet._aulPoints[2] = 2; - faces.push_back(aFacet); - //Fill the Kernel with the temp smesh structure and delete the current containers - aMesh.Adopt(vertices,faces); - - /////////////////////////////////////////////////////////////////////////// - //Now do Monte Carlo to minimize the BBox of the part - //Use Quaternions for the rotation stuff - Base::Rotation rotatex,rotatey,rotatez; - 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); - - //Rotate around each axes and choose the settings for the min bbox - Base::Matrix4D final_trafo; - Base::BoundBox3f aBBox,min_bbox; - double volumeBBOX,min_volumeBBOX; - //Get the current min_volumeBBOX and look if we find a lower one - aBBox = aMesh.GetBoundBox(); - min_volumeBBOX = aBBox.LengthX()*aBBox.LengthY()*aBBox.LengthZ(); - - MeshCore::MeshKernel atempkernel; - - float it_steps=10.0; - double step_size; - double alpha_x=0.0,alpha_y=0.0,alpha_z=0.0; - double perfect_ax=0.0,perfect_ay=0.0,perfect_az=0.0; - - //Do a Monte Carlo approach and start from the Principal Axis System - //and rotate +/- 60° around each axis in a first iteration - double angle_range_min_x=-M_PI/3.0,angle_range_max_x=M_PI/3.0, - angle_range_min_y=-M_PI/3.0,angle_range_max_y=M_PI/3.0, - angle_range_min_z=-M_PI/3.0,angle_range_max_z=M_PI/3.0; - - //We rotate until we are 0.1° sure to be in the right position - for (step_size = (2.0*M_PI/it_steps);step_size>(2.0*M_PI/3600.0);step_size=(2.0*M_PI/it_steps)) - { - for(alpha_x=angle_range_min_x;alpha_xgetFemMeshPtr()->transformGeometry(final_trafo); - ////////////////////////////////////////////////////////////////////////////////////// - //Now lets also do the movement to the 1st Quadrant in this function - aBBox = aMesh.GetBoundBox(); - //Get Distance vector from current lower left corner of BBox to origin - Base::Vector3f dist_vector; - dist_vector.x = -aBBox.MinX;dist_vector.y=-aBBox.MinY;dist_vector.z=-aBBox.MinZ; - Base::Matrix4D trans_matrix( - float(1.0),float(0.0),float(0.0),dist_vector.x, - float(0.0),float(1.0),float(0.0),dist_vector.y, - float(0.0),float(0.0),float(1.0),dist_vector.z, - float(0.0),float(0.0),float(0.0),float(1.0)); - inputMesh->getFemMeshPtr()->transformGeometry(trans_matrix); - - //inputMesh->getFemMeshPtr()->getSMesh()->ExportUNV("C:/temp/fine_tuning.unv"); - - } PY_CATCH; - - Py_Return; -} +//static PyObject * SMESH_PCA(PyObject *self, PyObject *args) +//{ +// PyObject *input; +// +// if (!PyArg_ParseTuple(args, "O",&input)) +// return NULL; +// +// PY_TRY { +// +// FemMeshPy *inputMesh = static_cast(input); +// MeshCore::MeshKernel aMesh; +// MeshCore::MeshPointArray vertices; +// vertices.clear(); +// MeshCore::MeshFacetArray faces; +// faces.clear(); +// MeshCore::MeshPoint current_node; +// SMDS_NodeIteratorPtr aNodeIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->nodesIterator(); +// for (;aNodeIter->more();) { +// const SMDS_MeshNode* aNode = aNodeIter->next(); +// current_node.Set(float(aNode->X()),float(aNode->Y()),float(aNode->Z())); +// vertices.push_back(current_node); +// } +// +// MeshCore::MeshFacet aFacet; +// aFacet._aulPoints[0] = 0;aFacet._aulPoints[1] = 1;aFacet._aulPoints[2] = 2; +// faces.push_back(aFacet); +// //Fill the Kernel with the temp smesh structure and delete the current containers +// aMesh.Adopt(vertices,faces); +// MeshCore::MeshEigensystem pca(aMesh); +// pca.Evaluate(); +// Base::Matrix4D Trafo = pca.Transform(); +// /*Let´s transform the input mesh with the PCA Matrix*/ +// inputMesh->getFemMeshPtr()->transformGeometry(Trafo); +// //inputMesh->getFemMeshPtr()->getSMesh()->ExportUNV("C:/Temp/PCA_alignment.unv"); +// //Now lets check if the smallest dimension of the BBox is oriented towards the Z-Axis. If not, lets rotate it around the X or Y axis +// //Use the SMESH structure for that +// // aMesh.Transform(Trafo); +// +// //Base::Rotation rotatex,rotatey,rotatez; +// //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); +// //double bbox_length_x,bbox_length_y,bbox_length_z; +// ////Rotate around the each axes and choose the settings for the min bbox +// //Base::Matrix4D final_trafo; +// //Base::BoundBox3f aBBox; +// ////Get the current BBOX and look for the size +// //aBBox = aMesh.GetBoundBox(); +// //bbox_length_x = aBBox.LengthX();bbox_length_y = aBBox.LengthY();bbox_length_z = aBBox.LengthZ(); +// ////Now do the rotation stuff +// //if (bbox_length_z < bbox_length_x && bbox_length_z < bbox_length_y) +// // Py_Return; +// //else if ( +// +// +// //MeshCore::MeshKernel atempkernel; +// +// //float it_steps=10.0; +// //double step_size; +// //double alpha_x=0.0,alpha_y=0.0,alpha_z=0.0; +// //double perfect_ax=0.0,perfect_ay=0.0,perfect_az=0.0; +// +// ////Do a Monte Carlo approach and start from the Principal Axis System +// ////and rotate +/- 60° around each axis in a first iteration +// //double angle_range_min_x=-PI/3.0,angle_range_max_x=PI/3.0, +// // angle_range_min_y=-PI/3.0,angle_range_max_y=PI/3.0, +// // angle_range_min_z=-PI/3.0,angle_range_max_z=PI/3.0; +// +// ////We rotate until we are 0.1° sure to be in the right position +// //for (step_size = (2.0*PI/it_steps);step_size>(2.0*PI/3600.0);step_size=(2.0*PI/it_steps)) +// //{ +// // for(alpha_x=angle_range_min_x;alpha_x(input); +// +// SMDS_VolumeIteratorPtr aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator(); +// Base::Vector3d a,b,c,a_b_product,temp,temp1; +// double volume =0.0; +// for (;aVolIter->more();) +// { +// const SMDS_MeshVolume* aVol = aVolIter->next(); +// //To make sure that the volume calculation is based on the ABAQUS element convention +// //The following Node mapping from SMESH to ABAQUS is necessary +// //ABAQUS_Node_Number|SMESH_Node_Number +// //0|0 +// //1|2 +// //2|1 +// //3|3 +// //4|6 +// //5|5 +// //6|4 +// //7|8 +// //8|9 +// //9|7 +// //The following coordinates of the little pyramids are based on ABAQUS convention and are numbered from +// //1 to 10 +// //1,5,8,7 +// temp.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); +// temp1.Set(aVol->GetNode(0)->X(),aVol->GetNode(0)->Y(),aVol->GetNode(0)->Z()); +// a = temp - temp1; +// temp.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); +// temp1.Set(aVol->GetNode(0)->X(),aVol->GetNode(0)->Y(),aVol->GetNode(0)->Z()); +// b = temp - temp1; +// temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); +// temp1.Set(aVol->GetNode(0)->X(),aVol->GetNode(0)->Y(),aVol->GetNode(0)->Z()); +// c = temp - temp1; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// //5,9,8,7 +// temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); +// temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); +// a = temp - temp1; +// temp.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); +// temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); +// b = temp - temp1; +// temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); +// temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); +// c = temp - temp1; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// //5,2,9,7 +// temp.Set(aVol->GetNode(2)->X(),aVol->GetNode(2)->Y(),aVol->GetNode(2)->Z()); +// temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); +// a = temp - temp1; +// temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); +// temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); +// b = temp - temp1; +// temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); +// temp1.Set(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); +// c = temp - temp1; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// //2,6,9,7 +// temp.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z()); +// temp1.Set(aVol->GetNode(2)->X(),aVol->GetNode(2)->Y(),aVol->GetNode(2)->Z()); +// a = temp - temp1; +// temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); +// temp1.Set(aVol->GetNode(2)->X(),aVol->GetNode(2)->Y(),aVol->GetNode(2)->Z()); +// b = temp - temp1; +// temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); +// temp1.Set(aVol->GetNode(2)->X(),aVol->GetNode(2)->Y(),aVol->GetNode(2)->Z()); +// c = temp - temp1; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// //9,6,10,7 +// temp.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z()); +// temp1.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); +// a = temp - temp1; +// temp.Set(aVol->GetNode(7)->X(),aVol->GetNode(7)->Y(),aVol->GetNode(7)->Z()); +// temp1.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); +// b = temp - temp1; +// temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); +// temp1.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); +// c = temp - temp1; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// //6,3,10,7 +// temp.Set(aVol->GetNode(1)->X(),aVol->GetNode(1)->Y(),aVol->GetNode(1)->Z()); +// temp1.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z()); +// a = temp - temp1; +// temp.Set(aVol->GetNode(7)->X(),aVol->GetNode(7)->Y(),aVol->GetNode(7)->Z()); +// temp1.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z()); +// b = temp - temp1; +// temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); +// temp1.Set(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z()); +// c = temp - temp1; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// //8,9,10,7 +// temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); +// temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); +// a = temp - temp1; +// temp.Set(aVol->GetNode(7)->X(),aVol->GetNode(7)->Y(),aVol->GetNode(7)->Z()); +// temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); +// b = temp - temp1; +// temp.Set(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); +// temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); +// c = temp - temp1; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// //8,9,10,4 +// temp.Set(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); +// temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); +// a = temp - temp1; +// temp.Set(aVol->GetNode(7)->X(),aVol->GetNode(7)->Y(),aVol->GetNode(7)->Z()); +// temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); +// b = temp - temp1; +// temp.Set(aVol->GetNode(3)->X(),aVol->GetNode(3)->Y(),aVol->GetNode(3)->Z()); +// temp1.Set(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); +// c = temp - temp1; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// } +// Py::Float py_volume(volume); +// return Py::new_reference_to(py_volume); +// +// } PY_CATCH; +// +// Py_Return; +//} + +//static PyObject * checkBB(PyObject *self, PyObject *args) +//{ +// PyObject *input; +// PyObject* plm=0; +// float billet_thickness; +// bool oversize = false; +// +// if (!PyArg_ParseTuple(args, "O|O!f", &input,&(Base::PlacementPy::Type),&plm,&billet_thickness)) +// return NULL; +// +// try { +// Base::Placement* placement = 0; +// if (plm) { +// placement = static_cast(plm)->getPlacementPtr(); +// +// } +// Base::Vector3d current_node; +// Base::Matrix4D matrix = placement->toMatrix(); +// FemMeshPy *inputMesh = static_cast(input); +// SMDS_NodeIteratorPtr aNodeIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->nodesIterator(); +// for (;aNodeIter->more();) { +// const SMDS_MeshNode* aNode = aNodeIter->next(); +// current_node.Set(float(aNode->X()),float(aNode->Y()),float(aNode->Z())); +// current_node = matrix * current_node; +// if(current_node.z > billet_thickness || current_node.z < -0.1) +// { +// //lets jump out of the function as soon as we find a +// //Node that is higher or lower than billet thickness +// oversize = true; +// Py::Boolean py_oversize(oversize); +// return Py::new_reference_to(py_oversize); +// } +// } +// Py::Boolean py_oversize(oversize); +// return Py::new_reference_to(py_oversize); +// } +// catch (const std::exception& e) { +// PyErr_SetString(PyExc_Exception, e.what()); +// return 0; +// } +// Py_Return; +//} +// +// +// +// +//static PyObject * getBoundary_Conditions(PyObject *self, PyObject *args) +//{ +// PyObject *input; +// Py::List boundary_nodes; +// +// if (!PyArg_ParseTuple(args, "O!", &(FemMeshPy::Type), &input)) +// //if (!PyArg_ParseTuple(args, "O",&input)) +// return NULL; +// +// PY_TRY { +// +// FemMeshPy *inputMesh = static_cast(input); +// MeshCore::MeshKernel aMesh; +// MeshCore::MeshPointArray vertices; +// vertices.clear(); +// MeshCore::MeshFacetArray faces; +// faces.clear(); +// MeshCore::MeshPoint current_node; +// SMDS_NodeIteratorPtr aNodeIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->nodesIterator(); +// SMDS_VolumeIteratorPtr aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator(); +// for (;aNodeIter->more();) { +// const SMDS_MeshNode* aNode = aNodeIter->next(); +// current_node.Set(float(aNode->X()),float(aNode->Y()),float(aNode->Z())); +// vertices.push_back(current_node); +// } +// MeshCore::MeshFacet aFacet; +// aFacet._aulPoints[0] = 0;aFacet._aulPoints[1] = 1;aFacet._aulPoints[2] = 2; +// faces.push_back(aFacet); +// //Fill the Kernel with the temp smesh structure and delete the current containers +// aMesh.Adopt(vertices,faces); +// Base::BoundBox3f aBBox; +// aBBox = aMesh.GetBoundBox(); +// +// float dist_length; +// Base::Vector3f dist; +// int minNodeID,maxNodeID,midNodeID; +// dist_length = FLOAT_MAX; +// +// aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator(); +// //We only search in non midside nodes (equals the first four nodes of each element) +// for (;aVolIter->more();) +// { +// const SMDS_MeshVolume* aVol = aVolIter->next(); +// for (unsigned j=0;j<4;j++) +// { +// const SMDS_MeshNode* aNode = aVol->GetNode(j); +// //Calc distance between the lower left corner and the most next point of the mesh +// dist.x = float(aNode->X())-aBBox.MinX;dist.y = float(aNode->Y())-aBBox.MinY;dist.z = float(aNode->Z())-aBBox.MinZ; +// if(dist.Length()GetID(); +// dist_length = dist.Length(); +// } +// } +// } +// +// boundary_nodes.append(Py::Int(minNodeID)); +// +// dist_length = FLOAT_MAX; +// aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator(); +// for (;aVolIter->more();) +// { +// const SMDS_MeshVolume* aVol = aVolIter->next(); +// for (unsigned j=0;j<4;j++) +// { +// const SMDS_MeshNode* aNode = aVol->GetNode(j); +// //Calc distance between the lower right corner and the most next point of the mesh +// dist.x = float(aNode->X())-aBBox.MaxX;dist.y = float(aNode->Y())-aBBox.MinY;dist.z = float(aNode->Z())-aBBox.MinZ; +// if(dist.Length()GetID(); +// dist_length = dist.Length(); +// } +// } +// } +// boundary_nodes.append(Py::Int(midNodeID)); +// +// dist_length = FLOAT_MAX; +// aVolIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->volumesIterator(); +// for (;aVolIter->more();) +// { +// const SMDS_MeshVolume* aVol = aVolIter->next(); +// for (unsigned j=0;j<4;j++) +// { +// const SMDS_MeshNode* aNode = aVol->GetNode(j); +// //Calc distance between the lowest lower right corner and the most next point of the mesh +// dist.x = float(aNode->X())-aBBox.MinX;dist.y = float(aNode->Y())-aBBox.MaxY;dist.z = float(aNode->Z())-aBBox.MinZ; +// if(dist.Length()GetID(); +// dist_length = dist.Length(); +// } +// } +// } +// boundary_nodes.append(Py::Int(maxNodeID)); +// +// +// +// +// return Py::new_reference_to(boundary_nodes); +// +// +// } PY_CATCH; +// +// Py_Return; +//} + + + + +//static PyObject * minBoundingBox(PyObject *self, PyObject *args) +//{ +// +// PyObject *input; +// +// if (!PyArg_ParseTuple(args, "O",&input)) +// return NULL; +// +// PY_TRY { +// FemMeshPy *inputMesh = static_cast(input); +// MeshCore::MeshKernel aMesh; +// MeshCore::MeshPointArray vertices; +// vertices.clear(); +// MeshCore::MeshFacetArray faces; +// faces.clear(); +// MeshCore::MeshPoint current_node; +// SMDS_NodeIteratorPtr aNodeIter = inputMesh->getFemMeshPtr()->getSMesh()->GetMeshDS()->nodesIterator(); +// for (;aNodeIter->more();) { +// const SMDS_MeshNode* aNode = aNodeIter->next(); +// current_node.Set(float(aNode->X()),float(aNode->Y()),float(aNode->Z())); +// vertices.push_back(current_node); +// } +// MeshCore::MeshFacet aFacet; +// aFacet._aulPoints[0] = 0;aFacet._aulPoints[1] = 1;aFacet._aulPoints[2] = 2; +// faces.push_back(aFacet); +// //Fill the Kernel with the temp smesh structure and delete the current containers +// aMesh.Adopt(vertices,faces); +// +// /////////////////////////////////////////////////////////////////////////// +// //Now do Monte Carlo to minimize the BBox of the part +// //Use Quaternions for the rotation stuff +// Base::Rotation rotatex,rotatey,rotatez; +// 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); +// +// //Rotate around each axes and choose the settings for the min bbox +// Base::Matrix4D final_trafo; +// Base::BoundBox3f aBBox,min_bbox; +// double volumeBBOX,min_volumeBBOX; +// //Get the current min_volumeBBOX and look if we find a lower one +// aBBox = aMesh.GetBoundBox(); +// min_volumeBBOX = aBBox.LengthX()*aBBox.LengthY()*aBBox.LengthZ(); +// +// MeshCore::MeshKernel atempkernel; +// +// float it_steps=10.0; +// double step_size; +// double alpha_x=0.0,alpha_y=0.0,alpha_z=0.0; +// double perfect_ax=0.0,perfect_ay=0.0,perfect_az=0.0; +// +// //Do a Monte Carlo approach and start from the Principal Axis System +// //and rotate +/- 60° around each axis in a first iteration +// double angle_range_min_x=-M_PI/3.0,angle_range_max_x=M_PI/3.0, +// angle_range_min_y=-M_PI/3.0,angle_range_max_y=M_PI/3.0, +// angle_range_min_z=-M_PI/3.0,angle_range_max_z=M_PI/3.0; +// +// //We rotate until we are 0.1° sure to be in the right position +// for (step_size = (2.0*M_PI/it_steps);step_size>(2.0*M_PI/3600.0);step_size=(2.0*M_PI/it_steps)) +// { +// for(alpha_x=angle_range_min_x;alpha_xgetFemMeshPtr()->transformGeometry(final_trafo); +// ////////////////////////////////////////////////////////////////////////////////////// +// //Now lets also do the movement to the 1st Quadrant in this function +// aBBox = aMesh.GetBoundBox(); +// //Get Distance vector from current lower left corner of BBox to origin +// Base::Vector3f dist_vector; +// dist_vector.x = -aBBox.MinX;dist_vector.y=-aBBox.MinY;dist_vector.z=-aBBox.MinZ; +// Base::Matrix4D trans_matrix( +// float(1.0),float(0.0),float(0.0),dist_vector.x, +// float(0.0),float(1.0),float(0.0),dist_vector.y, +// float(0.0),float(0.0),float(1.0),dist_vector.z, +// float(0.0),float(0.0),float(0.0),float(1.0)); +// inputMesh->getFemMeshPtr()->transformGeometry(trans_matrix); +// +// //inputMesh->getFemMeshPtr()->getSMesh()->ExportUNV("C:/temp/fine_tuning.unv"); +// +// } PY_CATCH; +// +// Py_Return; +//} static PyObject * importer(PyObject *self, PyObject *args) @@ -705,360 +705,360 @@ static PyObject * importer(PyObject *self, PyObject *args) } -static PyObject * import_NASTRAN(PyObject *self, PyObject *args) -{ - const char* filename_input, *filename_output; - if (!PyArg_ParseTuple(args, "ss",&filename_input,&filename_output)) - return NULL; - - PY_TRY { - - std::ifstream inputfile; - - //Für Debugoutput - ofstream anOutput; - anOutput.open("c:/time_measurement.txt"); - time_t seconds1,seconds2; - - inputfile.open(filename_input); - if (!inputfile.is_open()) //Exists...? - { - std::cerr << "File not found. Exiting..." << std::endl; - return NULL; - } - - //Return the line pointer to the beginning of the file - inputfile.seekg(std::ifstream::beg); - std::string line1,line2,temp; - std::stringstream astream; - std::vector tetra_element; - std::vector element_id; - element_id.clear(); - std::vector > all_elements; - std::vector >::iterator all_element_iterator; - std::vector::iterator one_element_iterator; - all_elements.clear(); - MeshCore::MeshKernel aMesh; - MeshCore::MeshPointArray vertices; - vertices.clear(); - MeshCore::MeshFacetArray faces; - faces.clear(); - MeshCore::MeshPoint current_node; - - seconds1 = time(NULL); - do - { - std::getline(inputfile,line1); - if (line1.size() == 0) continue; - if (line1.find("GRID*")!= std::string::npos) //We found a Grid line - { - //Now lets extract the GRID Points = Nodes - //As each GRID Line consists of two subsequent lines we have to - //take care of that as well - std::getline(inputfile,line2); - //Extract X Value - current_node.x = float(atof(line1.substr(40,56).c_str())); - //Extract Y Value - current_node.y = float(atof(line1.substr(56,72).c_str())); - //Extract Z Value - current_node.z = float(atof(line2.substr(8,24).c_str())); - - vertices.push_back(current_node); - } - else if (line1.find("CTETRA")!= std::string::npos) - { - tetra_element.clear(); - //Lets extract the elements - //As each Element Line consists of two subsequent lines as well - //we have to take care of that - //At a first step we only extract Quadratic Tetrahedral Elements - std::getline(inputfile,line2); - element_id.push_back(atoi(line1.substr(8,16).c_str())); - tetra_element.push_back(atoi(line1.substr(24,32).c_str())); - tetra_element.push_back(atoi(line1.substr(32,40).c_str())); - tetra_element.push_back(atoi(line1.substr(40,48).c_str())); - tetra_element.push_back(atoi(line1.substr(48,56).c_str())); - tetra_element.push_back(atoi(line1.substr(56,64).c_str())); - tetra_element.push_back(atoi(line1.substr(64,72).c_str())); - tetra_element.push_back(atoi(line2.substr(8,16).c_str())); - tetra_element.push_back(atoi(line2.substr(16,24).c_str())); - tetra_element.push_back(atoi(line2.substr(24,32).c_str())); - tetra_element.push_back(atoi(line2.substr(32,40).c_str())); - - all_elements.push_back(tetra_element); - } - - } - while (inputfile.good()); - inputfile.close(); - - seconds2 = time(NULL); - - anOutput << seconds2-seconds1 <<" for Parsing the input file"<(2.0*M_PI/360.0);step_size=(2.0*M_PI/it_steps)) - { - for(alpha_x=angle_range_min_x;alpha_xCreateMesh(1, true); - SMESHDS_Mesh* meshds = mesh->GetMeshDS(); - - MeshCore::MeshPointIterator anodeiterator(aMesh); - int j=1; - for(anodeiterator.Begin(); anodeiterator.More(); anodeiterator.Next()) - { - meshds->AddNodeWithID((*anodeiterator).x,(*anodeiterator).y,(*anodeiterator).z,j); - j++; - } - - for(unsigned int i=0;iAddVolumeWithID( - meshds->FindNode(all_elements[i][0]), - meshds->FindNode(all_elements[i][2]), - meshds->FindNode(all_elements[i][1]), - meshds->FindNode(all_elements[i][3]), - meshds->FindNode(all_elements[i][6]), - meshds->FindNode(all_elements[i][5]), - meshds->FindNode(all_elements[i][4]), - meshds->FindNode(all_elements[i][9]), - meshds->FindNode(all_elements[i][7]), - meshds->FindNode(all_elements[i][8]), - element_id[i] - ); - } - - mesh->ExportUNV(filename_output); - ////////////////////////////////////////////////////////////////////////////////////////// - seconds2 = time(NULL); - anOutput << seconds2-seconds1 << " seconds for the Mesh Export" << std::endl; - - //Output also to ABAQUS Input Format - ofstream anABAQUS_Output; - anABAQUS_Output.open("d:/abaqus_output.inp"); - anABAQUS_Output << "*Node , NSET=Nall" << std::endl; - j=1; - for(anodeiterator.Begin(); anodeiterator.More(); anodeiterator.Next()) - { - anABAQUS_Output << j <<"," - <<(*anodeiterator).x << "," - <<(*anodeiterator).y << "," - <<(*anodeiterator).z << std::endl; - j++; - } - anABAQUS_Output << "*Element, TYPE=C3D10, ELSET=Eall" << std::endl; - j=1; - for(unsigned int i=0;iat(4)-1]-vertices[all_element_iterator->at(0)-1]; - b = vertices[all_element_iterator->at(7)-1]-vertices[all_element_iterator->at(0)-1]; - c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(0)-1]; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - //5,9,8,7 - a = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(4)-1]; - b = vertices[all_element_iterator->at(7)-1]-vertices[all_element_iterator->at(4)-1]; - c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(4)-1]; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - //5,2,9,7 - a = vertices[all_element_iterator->at(1)-1]-vertices[all_element_iterator->at(4)-1]; - b = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(4)-1]; - c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(4)-1]; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - //2,6,9,7 - a = vertices[all_element_iterator->at(5)-1]-vertices[all_element_iterator->at(1)-1]; - b = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(1)-1]; - c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(1)-1]; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - //9,6,10,7 - a = vertices[all_element_iterator->at(5)-1]-vertices[all_element_iterator->at(8)-1]; - b = vertices[all_element_iterator->at(9)-1]-vertices[all_element_iterator->at(8)-1]; - c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(8)-1]; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - //6,3,10,7 - a = vertices[all_element_iterator->at(2)-1]-vertices[all_element_iterator->at(5)-1]; - b = vertices[all_element_iterator->at(9)-1]-vertices[all_element_iterator->at(5)-1]; - c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(5)-1]; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - //8,9,10,7 - a = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(7)-1]; - b = vertices[all_element_iterator->at(9)-1]-vertices[all_element_iterator->at(7)-1]; - c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(7)-1]; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - //8,9,10,4 - a = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(7)-1]; - b = vertices[all_element_iterator->at(9)-1]-vertices[all_element_iterator->at(7)-1]; - c = vertices[all_element_iterator->at(3)-1]-vertices[all_element_iterator->at(7)-1]; - 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; - volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); - - } - seconds2=time(NULL); - anOutput << seconds2-seconds1 << " seconds for Volume Calculation " << "Volumen " << volume/1000.0/1000.0/1000.0 << " in m^3" << std::endl; - anOutput << "Volumen der BBox" << min_volumeBBOX/1000.0/1000.0/1000.0 << std::endl; - anOutput << "Fly to Buy Ratio: " << min_volumeBBOX / volume << std::endl; - anOutput.close(); - - } PY_CATCH; - - Py_Return; -} +//static PyObject * import_NASTRAN(PyObject *self, PyObject *args) +//{ +// const char* filename_input, *filename_output; +// if (!PyArg_ParseTuple(args, "ss",&filename_input,&filename_output)) +// return NULL; +// +// PY_TRY { +// +// std::ifstream inputfile; +// +// //Für Debugoutput +// ofstream anOutput; +// anOutput.open("c:/time_measurement.txt"); +// time_t seconds1,seconds2; +// +// inputfile.open(filename_input); +// if (!inputfile.is_open()) //Exists...? +// { +// std::cerr << "File not found. Exiting..." << std::endl; +// return NULL; +// } +// +// //Return the line pointer to the beginning of the file +// inputfile.seekg(std::ifstream::beg); +// std::string line1,line2,temp; +// std::stringstream astream; +// std::vector tetra_element; +// std::vector element_id; +// element_id.clear(); +// std::vector > all_elements; +// std::vector >::iterator all_element_iterator; +// std::vector::iterator one_element_iterator; +// all_elements.clear(); +// MeshCore::MeshKernel aMesh; +// MeshCore::MeshPointArray vertices; +// vertices.clear(); +// MeshCore::MeshFacetArray faces; +// faces.clear(); +// MeshCore::MeshPoint current_node; +// +// seconds1 = time(NULL); +// do +// { +// std::getline(inputfile,line1); +// if (line1.size() == 0) continue; +// if (line1.find("GRID*")!= std::string::npos) //We found a Grid line +// { +// //Now lets extract the GRID Points = Nodes +// //As each GRID Line consists of two subsequent lines we have to +// //take care of that as well +// std::getline(inputfile,line2); +// //Extract X Value +// current_node.x = float(atof(line1.substr(40,56).c_str())); +// //Extract Y Value +// current_node.y = float(atof(line1.substr(56,72).c_str())); +// //Extract Z Value +// current_node.z = float(atof(line2.substr(8,24).c_str())); +// +// vertices.push_back(current_node); +// } +// else if (line1.find("CTETRA")!= std::string::npos) +// { +// tetra_element.clear(); +// //Lets extract the elements +// //As each Element Line consists of two subsequent lines as well +// //we have to take care of that +// //At a first step we only extract Quadratic Tetrahedral Elements +// std::getline(inputfile,line2); +// element_id.push_back(atoi(line1.substr(8,16).c_str())); +// tetra_element.push_back(atoi(line1.substr(24,32).c_str())); +// tetra_element.push_back(atoi(line1.substr(32,40).c_str())); +// tetra_element.push_back(atoi(line1.substr(40,48).c_str())); +// tetra_element.push_back(atoi(line1.substr(48,56).c_str())); +// tetra_element.push_back(atoi(line1.substr(56,64).c_str())); +// tetra_element.push_back(atoi(line1.substr(64,72).c_str())); +// tetra_element.push_back(atoi(line2.substr(8,16).c_str())); +// tetra_element.push_back(atoi(line2.substr(16,24).c_str())); +// tetra_element.push_back(atoi(line2.substr(24,32).c_str())); +// tetra_element.push_back(atoi(line2.substr(32,40).c_str())); +// +// all_elements.push_back(tetra_element); +// } +// +// } +// while (inputfile.good()); +// inputfile.close(); +// +// seconds2 = time(NULL); +// +// anOutput << seconds2-seconds1 <<" for Parsing the input file"<(2.0*M_PI/360.0);step_size=(2.0*M_PI/it_steps)) +// { +// for(alpha_x=angle_range_min_x;alpha_xCreateMesh(1, true); +// SMESHDS_Mesh* meshds = mesh->GetMeshDS(); +// +// MeshCore::MeshPointIterator anodeiterator(aMesh); +// int j=1; +// for(anodeiterator.Begin(); anodeiterator.More(); anodeiterator.Next()) +// { +// meshds->AddNodeWithID((*anodeiterator).x,(*anodeiterator).y,(*anodeiterator).z,j); +// j++; +// } +// +// for(unsigned int i=0;iAddVolumeWithID( +// meshds->FindNode(all_elements[i][0]), +// meshds->FindNode(all_elements[i][2]), +// meshds->FindNode(all_elements[i][1]), +// meshds->FindNode(all_elements[i][3]), +// meshds->FindNode(all_elements[i][6]), +// meshds->FindNode(all_elements[i][5]), +// meshds->FindNode(all_elements[i][4]), +// meshds->FindNode(all_elements[i][9]), +// meshds->FindNode(all_elements[i][7]), +// meshds->FindNode(all_elements[i][8]), +// element_id[i] +// ); +// } +// +// mesh->ExportUNV(filename_output); +// ////////////////////////////////////////////////////////////////////////////////////////// +// seconds2 = time(NULL); +// anOutput << seconds2-seconds1 << " seconds for the Mesh Export" << std::endl; +// +// //Output also to ABAQUS Input Format +// ofstream anABAQUS_Output; +// anABAQUS_Output.open("d:/abaqus_output.inp"); +// anABAQUS_Output << "*Node , NSET=Nall" << std::endl; +// j=1; +// for(anodeiterator.Begin(); anodeiterator.More(); anodeiterator.Next()) +// { +// anABAQUS_Output << j <<"," +// <<(*anodeiterator).x << "," +// <<(*anodeiterator).y << "," +// <<(*anodeiterator).z << std::endl; +// j++; +// } +// anABAQUS_Output << "*Element, TYPE=C3D10, ELSET=Eall" << std::endl; +// j=1; +// for(unsigned int i=0;iat(4)-1]-vertices[all_element_iterator->at(0)-1]; +// b = vertices[all_element_iterator->at(7)-1]-vertices[all_element_iterator->at(0)-1]; +// c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(0)-1]; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// //5,9,8,7 +// a = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(4)-1]; +// b = vertices[all_element_iterator->at(7)-1]-vertices[all_element_iterator->at(4)-1]; +// c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(4)-1]; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// //5,2,9,7 +// a = vertices[all_element_iterator->at(1)-1]-vertices[all_element_iterator->at(4)-1]; +// b = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(4)-1]; +// c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(4)-1]; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// //2,6,9,7 +// a = vertices[all_element_iterator->at(5)-1]-vertices[all_element_iterator->at(1)-1]; +// b = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(1)-1]; +// c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(1)-1]; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// //9,6,10,7 +// a = vertices[all_element_iterator->at(5)-1]-vertices[all_element_iterator->at(8)-1]; +// b = vertices[all_element_iterator->at(9)-1]-vertices[all_element_iterator->at(8)-1]; +// c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(8)-1]; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// //6,3,10,7 +// a = vertices[all_element_iterator->at(2)-1]-vertices[all_element_iterator->at(5)-1]; +// b = vertices[all_element_iterator->at(9)-1]-vertices[all_element_iterator->at(5)-1]; +// c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(5)-1]; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// //8,9,10,7 +// a = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(7)-1]; +// b = vertices[all_element_iterator->at(9)-1]-vertices[all_element_iterator->at(7)-1]; +// c = vertices[all_element_iterator->at(6)-1]-vertices[all_element_iterator->at(7)-1]; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// //8,9,10,4 +// a = vertices[all_element_iterator->at(8)-1]-vertices[all_element_iterator->at(7)-1]; +// b = vertices[all_element_iterator->at(9)-1]-vertices[all_element_iterator->at(7)-1]; +// c = vertices[all_element_iterator->at(3)-1]-vertices[all_element_iterator->at(7)-1]; +// 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; +// volume += 1.0/6.0 * fabs((a_b_product.x * c.x)+ (a_b_product.y * c.y)+(a_b_product.z * c.z)); +// +// } +// seconds2=time(NULL); +// anOutput << seconds2-seconds1 << " seconds for Volume Calculation " << "Volumen " << volume/1000.0/1000.0/1000.0 << " in m^3" << std::endl; +// anOutput << "Volumen der BBox" << min_volumeBBOX/1000.0/1000.0/1000.0 << std::endl; +// anOutput << "Fly to Buy Ratio: " << min_volumeBBOX / volume << std::endl; +// anOutput.close(); +// +// } PY_CATCH; +// +// Py_Return; +//} static PyObject * exporter(PyObject *self, PyObject *args) @@ -1105,14 +1105,14 @@ struct PyMethodDef Fem_methods[] = { {"read" ,read, Py_NEWARGS, "Read a mesh from a file and returns a Mesh object."}, {"show" ,show ,METH_VARARGS, "show(shape) -- Add the shape to the active document or create one if no document exists."}, - {"calcMeshVolume", calcMeshVolume, Py_NEWARGS, - "Calculate Mesh Volume for C3D10"}, - {"getBoundary_Conditions" , getBoundary_Conditions, Py_NEWARGS, - "Get Boundary Conditions for Residual Stress Calculation"}, - {"SMESH_PCA" , SMESH_PCA, Py_NEWARGS, - "Get a Matrix4D related to the PCA of a Mesh Object"}, - {"import_NASTRAN",import_NASTRAN, Py_NEWARGS, "Import Nastran files, for tests only. Use read() or insert()"}, - {"minBoundingBox",minBoundingBox,Py_NEWARGS,"Minimize the Bounding Box and reorient the mesh to the 1st Quadrant"}, - {"checkBB",checkBB,Py_NEWARGS,"Check if the nodal z-values are still in the prescribed range"}, + // {"calcMeshVolume", calcMeshVolume, Py_NEWARGS, + // "Calculate Mesh Volume for C3D10"}, + // {"getBoundary_Conditions" , getBoundary_Conditions, Py_NEWARGS, + // "Get Boundary Conditions for Residual Stress Calculation"}, + //{"SMESH_PCA" , SMESH_PCA, Py_NEWARGS, + // "Get a Matrix4D related to the PCA of a Mesh Object"}, + //{"import_NASTRAN",import_NASTRAN, Py_NEWARGS, "Import Nastran files, for tests only. Use read() or insert()"}, + //{"minBoundingBox",minBoundingBox,Py_NEWARGS,"Minimize the Bounding Box and reorient the mesh to the 1st Quadrant"}, + //{"checkBB",checkBB,Py_NEWARGS,"Check if the nodal z-values are still in the prescribed range"}, {NULL, NULL} /* sentinel */ }; diff --git a/src/Mod/Fem/App/CMakeLists.txt b/src/Mod/Fem/App/CMakeLists.txt index f4cc6dd98..944cca71d 100755 --- a/src/Mod/Fem/App/CMakeLists.txt +++ b/src/Mod/Fem/App/CMakeLists.txt @@ -28,7 +28,6 @@ link_directories(${OCC_LIBRARY_DIR}) if(FREECAD_BUILD_FEM_NETGEN) set(Fem_LIBS Part - Mesh FreeCADApp StdMeshers NETGENPlugin @@ -37,7 +36,6 @@ if(FREECAD_BUILD_FEM_NETGEN) else(FREECAD_BUILD_FEM_NETGEN) set(Fem_LIBS Part - Mesh FreeCADApp StdMeshers SMESH