Volum for Tet10 FemMesh

FreeCAD · Dec 9, 2013 · 132b75e · 132b75e
1 parent 401cb13
commit 132b75e
Show file tree

Hide file tree

Showing 4 changed files with 123 additions and 6 deletions.
diff --git 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;i<all_elements.size();i++)
+//		{
+//			//Die Reihenfolge wie hier die Elemente hinzugefügt werden ist sehr wichtig. 
+//			//Ansonsten ist eine konsistente Datenstruktur nicht möglich
+//			meshds->AddVolumeWithID(
+//				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);
+
+
+}
diff --git a/src/Mod/Fem/App/FemMesh.h b/src/Mod/Fem/App/FemMesh.h
@@ -26,6 +26,7 @@
 
 #include <App/ComplexGeoData.h>
 #include <Base/Placement.h>
+#include <Base/Quantity.h>
 
 #include <vector>
 #include <list>
@@ -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
@@ -179,11 +179,17 @@
 		  </Documentation>
 		  <Parameter Name="SubMeshCount" Type="Int"/>
 	  </Attribute>
-	  <Attribute Name="GroupCount" ReadOnly="true">
-		  <Documentation>
-			  <UserDocu>Number of Groups in the Mesh.</UserDocu>
-		  </Documentation>
-		  <Parameter Name="GroupCount" Type="Int"/>
-	  </Attribute>
+      <Attribute Name="GroupCount" ReadOnly="true">
+          <Documentation>
+              <UserDocu>Number of Groups in the Mesh.</UserDocu>
+          </Documentation>
+          <Parameter Name="GroupCount" Type="Int"/>
+      </Attribute>
+      <Attribute Name="Volume" ReadOnly="true">
+          <Documentation>
+              <UserDocu>Volume of the mesh.</UserDocu>
+          </Documentation>
+          <Parameter Name="Volume" Type="Object"/>
+      </Attribute>
   </PythonExport>
 </GenerateModel>
diff --git a/src/Mod/Fem/App/FemMeshPyImp.cpp b/src/Mod/Fem/App/FemMeshPyImp.cpp
@@ -34,6 +34,7 @@
 #include <Base/VectorPy.h>
 #include <Base/MatrixPy.h>
 #include <Base/PlacementPy.h>
+#include <Base/QuantityPy.h>
 
 #include <Mod/Part/App/TopoShapePy.h>
 #include <Mod/Part/App/TopoShapeFacePy.h>
@@ -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