From 132b75e5ab714710dd56ae3b95e489f7bf66cfa5 Mon Sep 17 00:00:00 2001 From: jriegel Date: Mon, 9 Dec 2013 23:48:49 +0100 Subject: [PATCH] Volum for Tet10 FemMesh --- src/Mod/Fem/App/FemMesh.cpp | 102 +++++++++++++++++++++++++++++++ src/Mod/Fem/App/FemMesh.h | 3 + src/Mod/Fem/App/FemMeshPy.xml | 18 ++++-- src/Mod/Fem/App/FemMeshPyImp.cpp | 6 ++ 4 files changed, 123 insertions(+), 6 deletions(-) diff --git a/src/Mod/Fem/App/FemMesh.cpp b/src/Mod/Fem/App/FemMesh.cpp index 0d3323549e61..02eebd515802 100755 --- a/src/Mod/Fem/App/FemMesh.cpp +++ b/src/Mod/Fem/App/FemMesh.cpp @@ -891,3 +891,105 @@ struct Fem::FemMesh::FemMeshInfo FemMesh::getInfo(void) const{ return rtrn; } +// 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] +// ); +// } + +Base::Quantity FemMesh::getVolume(void)const +{ + SMDS_VolumeIteratorPtr aVolIter = myMesh->GetMeshDS()->volumesIterator(); + + //Calculate Mesh Volume + //For an accurate Volume Calculation of a quadratic Tetrahedron + //we have to calculate the Volume of 8 Sub-Tetrahedrons + Base::Vector3d a,b,c,a_b_product; + double volume = 0.0; + + for (;aVolIter->more();) + { + const SMDS_MeshVolume* aVol = aVolIter->next(); + + if ( aVol->NbNodes() != 10 ) continue; + + Base::Vector3d v1(aVol->GetNode(1)->X(),aVol->GetNode(1)->Y(),aVol->GetNode(1)->Z()); + Base::Vector3d v0(aVol->GetNode(0)->X(),aVol->GetNode(0)->Y(),aVol->GetNode(0)->Z()); + Base::Vector3d v2(aVol->GetNode(2)->X(),aVol->GetNode(2)->Y(),aVol->GetNode(2)->Z()); + Base::Vector3d v3(aVol->GetNode(3)->X(),aVol->GetNode(3)->Y(),aVol->GetNode(3)->Z()); + Base::Vector3d v4(aVol->GetNode(4)->X(),aVol->GetNode(4)->Y(),aVol->GetNode(4)->Z()); + Base::Vector3d v6(aVol->GetNode(6)->X(),aVol->GetNode(6)->Y(),aVol->GetNode(6)->Z()); + Base::Vector3d v5(aVol->GetNode(5)->X(),aVol->GetNode(5)->Y(),aVol->GetNode(5)->Z()); + Base::Vector3d v8(aVol->GetNode(8)->X(),aVol->GetNode(8)->Y(),aVol->GetNode(8)->Z()); + Base::Vector3d v7(aVol->GetNode(7)->X(),aVol->GetNode(7)->Y(),aVol->GetNode(7)->Z()); + Base::Vector3d v9(aVol->GetNode(9)->X(),aVol->GetNode(9)->Y(),aVol->GetNode(9)->Z()); + + + //1,5,8,7 + a = v4 -v0 ; + b = v7 -v0 ; + c = v6 -v0 ; + 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 = v8 -v4 ; + b = v7 -v4 ; + c = v6 -v4 ; + 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 = v1 -v4 ; + b = v8 -v4 ; + c = v6 -v4 ; + 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 = v5 -v1 ; + b = v8 -v1 ; + c = v6 -v1 ; + 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 = v5 -v8 ; + b = v9 -v8 ; + c = v6 -v8 ; + 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 = v2 -v5 ; + b = v9 -v5 ; + c = v6 -v5 ; + 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 = v8 -v7 ; + b = v9 -v7 ; + c = v6 -v7 ; + 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 = v8 -v7 ; + b = v9 -v7 ; + c = v3 -v7 ; + 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)); + + } + + return Base::Quantity(volume,Unit::Volume); + + +} \ No newline at end of file diff --git a/src/Mod/Fem/App/FemMesh.h b/src/Mod/Fem/App/FemMesh.h index 70c330636412..0d2cf3bed3fe 100755 --- a/src/Mod/Fem/App/FemMesh.h +++ b/src/Mod/Fem/App/FemMesh.h @@ -26,6 +26,7 @@ #include #include +#include #include #include @@ -96,6 +97,8 @@ class AppFemExport FemMesh : public Data::ComplexGeoData Base::Matrix4D getTransform(void) const; /// Bound box from the shape Base::BoundBox3d getBoundBox(void)const; + /// get the volume (when there are volume elements) + Base::Quantity getVolume(void)const; //@} /** @name Modification */ diff --git a/src/Mod/Fem/App/FemMeshPy.xml b/src/Mod/Fem/App/FemMeshPy.xml index ebc7b0b8fc9f..8933354f8865 100755 --- a/src/Mod/Fem/App/FemMeshPy.xml +++ b/src/Mod/Fem/App/FemMeshPy.xml @@ -179,11 +179,17 @@ - - - Number of Groups in the Mesh. - - - + + + Number of Groups in the Mesh. + + + + + + Volume of the mesh. + + + diff --git a/src/Mod/Fem/App/FemMeshPyImp.cpp b/src/Mod/Fem/App/FemMeshPyImp.cpp index dcae203bb5a0..0cfef63ee0fb 100755 --- a/src/Mod/Fem/App/FemMeshPyImp.cpp +++ b/src/Mod/Fem/App/FemMeshPyImp.cpp @@ -34,6 +34,7 @@ #include #include #include +#include #include #include @@ -603,6 +604,11 @@ Py::Int FemMeshPy::getGroupCount(void) const return Py::Int(getFemMeshPtr()->getSMesh()->NbGroup()); } +Py::Object FemMeshPy::getVolume(void) const +{ + return Py::Object(new Base::QuantityPy(new Base::Quantity(getFemMeshPtr()->getVolume()))); + +} // ===== custom attributes ============================================================ PyObject *FemMeshPy::getCustomAttributes(const char* /*attr*/) const