Area

src/GaussQuadrature.cpp

@@ -218,3 +218,63 @@ void GaussQuadrature::GetPoints(
 
 ///////////////////////////////////////////////////////////////////////////////
 
+void GaussQuadrature::GetPoints(
+	int nCount,
+	DataArray2D<double> & dG,
+	DataArray1D<double> & dW
+) {
+	if(nCount <= 4){
+
+		dW[0] = 0.0549758718276609338191631624501052;
+        	dW[1] = 0.0549758718276609338191631624501052;
+        	dW[2] = 0.0549758718276609338191631624501052;
+        	dW[3] = 0.111690794839005732847503504216561;
+        	dW[4] = 0.111690794839005732847503504216561;
+        	dW[5] = 0.111690794839005732847503504216561;
+
+		dG[0][0] = 0.8168475729804585130808570731956, dG[0][1] = 0.0915762135097707434595714634022015;
+        	dG[1][0] = 0.0915762135097707434595714634022015, dG[1][1] = 0.8168475729804585130808570731956;
+        	dG[2][0] = 0.0915762135097707434595714634022015, dG[2][1] = 0.0915762135097707434595714634022015;
+        	dG[3][0] = 0.1081030181680702273633414922339, dG[3][1] = 0.445948490915964886318329253883051;
+        	dG[4][0] = 0.445948490915964886318329253883051, dG[4][1] = 0.1081030181680702273633414922339;
+        	dG[5][0] = 0.445948490915964886318329253883051, dG[5][1] = 0.445948490915964886318329253883051;
+	}
+	else if(nCount <= 8){
+		dW[0] = 0.072157803838893584125545555249701;
+        	dW[1] = 0.051608685267359125140895775145648;
+		dW[2] = 0.051608685267359125140895775145648;
+		dW[3] = 0.051608685267359125140895775145648;
+		dW[4] = 0.016229248811599040155462964170437;
+		dW[5] = 0.016229248811599040155462964170437;
+		dW[6] = 0.016229248811599040155462964170437;
+		dW[7] = 0.047545817133642312396948052190887;
+		dW[8] = 0.047545817133642312396948052190887;
+		dW[9] = 0.047545817133642312396948052190887;
+		dW[10] = 0.013615157087217497132422345038231;
+		dW[11] = 0.013615157087217497132422345038231;
+		dW[12] = 0.013615157087217497132422345038231;
+		dW[13] = 0.013615157087217497132422345038231;
+		dW[14] = 0.013615157087217497132422345038231;
+		dW[15] = 0.01361515708721749713242234503823;
+
+		dG[0][0] = 0.33333333333333333333333333333333, dG[0][1] = 0.33333333333333333333333333333333;
+		dG[1][0] = 0.1705693077517602066222935014994, dG[1][1] = 0.1705693077517602066222935014994;
+		dG[2][0] = 0.1705693077517602066222935014994, dG[2][1] = 0.65886138449647958675541299700121;
+		dG[3][0] = 0.65886138449647958675541299700121, dG[3][1] = 0.1705693077517602066222935014994;
+		dG[4][0] = 0.050547228317030975458423550596387, dG[4][1] = 0.050547228317030975458423550596387;
+		dG[5][0] = 0.050547228317030975458423550596387, dG[5][1] = 0.89890554336593804908315289880723;
+		dG[6][0] = 0.89890554336593804908315289880723, dG[6][1] = 0.050547228317030975458423550596387;
+		dG[7][0] = 0.45929258829272315602881551450124, dG[7][1] = 0.45929258829272315602881551450124;
+		dG[8][0] = 0.45929258829272315602881551450124, dG[8][1] = 0.081414823414553687942368970997513;
+		dG[9][0] = 0.081414823414553687942368970997513, dG[9][1] = 0.45929258829272315602881551450124;
+		dG[10][0] = 0.72849239295540428124100037918962, dG[10][1] = 0.26311282963463811342178578626121;
+		dG[11][0] = 0.26311282963463811342178578626121, dG[11][1] = 0.72849239295540428124100037918962;
+		dG[12][0] = 0.72849239295540428124100037918962, dG[12][1] = 0.0083947774099576053372138345491687;
+		dG[13][0] = 0.0083947774099576053372138345491687, dG[13][1] = 0.72849239295540428124100037918962;
+		dG[14][0] = 0.26311282963463811342178578626121, dG[14][1] = 0.0083947774099576053372138345491687;
+		dG[15][0] = 0.0083947774099576053372138345491687, dG[15][1] = 0.26311282963463811342178578626121;
+
+	}	
+
+}
+

src/GaussQuadrature.h

@@ -18,6 +18,7 @@
 #define _GAUSSQUADRATURE_H_
 
 #include "DataArray1D.h"
+#include "DataArray2D.h"
 
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -48,6 +49,17 @@ class GaussQuadrature {
 		DataArray1D<double> & dG,
 		DataArray1D<double> & dW
 	);
+
+	///	<summary>
+	///		Return the two-dimensional Gauss quadrature points and their corresponding
+	///		weights for the given number of points and reference element.
+	///	</summary>
+
+	static void GetPoints(
+		int nCount,
+		DataArray2D<double> & dG,
+		DataArray1D<double> & dW
+	);
 };
 
 ///////////////////////////////////////////////////////////////////////////////

src/GridElements.cpp

@@ -1895,6 +1895,194 @@ Real CalculateFaceAreaQuadratureMethod(
 
 ///////////////////////////////////////////////////////////////////////////////
 
+Real CalculateFaceAreaQuadratureMethod2(
+	const Face & face,
+	const NodeVector & nodes
+) {
+
+	int nOrder1 = 4;
+	int nOrder2 = 8;
+
+	double dFaceArea = 0.0;
+
+	double h = MaxEdgeLength(face,nodes);
+
+	if(h < 0.004){
+
+		dFaceArea = CalculateFaceAreaQuadrature(face,nodes,nOrder1);
+
+	}
+	else if(h < 0.09){
+
+		dFaceArea = CalculateFaceAreaQuadrature(face,nodes,nOrder2);
+
+	}
+	else{
+
+		FaceVector faces;
+		faces.push_back(face);
+		dFaceArea = CalculateFaceAreaQuadratureSplit(faces,nodes,nOrder2);
+
+	}
+
+	return dFaceArea;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+Real CalculateFaceAreaQuadrature(
+	const Face & face,
+	const NodeVector & nodes,
+	int nOrder
+) {
+
+	double dFaceArea=0.0;
+	DataArray2D<double> dG1;
+	DataArray1D<double> dW1;
+	GaussQuadrature::GetPoints(nOrder, dG1, dW1);
+
+	int nTriangles = face.edges.size() - 2;
+
+	for (int j = 0; j < nTriangles; j++) {
+
+		Node node1 = nodes[face[0]];
+		Node node2 = nodes[face[j+1]];
+		Node node3 = nodes[face[j+2]];
+
+		Node nodeCross = CrossProduct(node2,node3);
+
+		double tri_pro = node1.x*nodeCross.x + node1.y*nodeCross.y + node1.z*nodeCross.z;
+
+		for (int q =0; q < dW1.GetRows(); q++){
+
+			Node pnts_q = node1*(1.0-dG1[q][0]-dG1[q][1]) + node2*dG1[q][0] + node3*dG1[q][1];
+
+			double nrm_q = pnts_q.Magnitude();
+
+			Node sph_q = pnts_q/nrm_q;
+
+			double nrm_q3 = pow(nrm_q,3);
+
+			dFaceArea += (tri_pro/nrm_q3)*dW1[q];
+		}
+
+	}
+
+	return dFaceArea;
+
+
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+Real CalculateFaceAreaQuadratureSplit(
+	const FaceVector & faces,
+	const NodeVector & nodes,
+	int & nOrder
+) {
+
+	double tol = 0.0;
+	double dFaceArea = 0.0;
+
+	if(nOrder >= 8){
+        	tol = 0.05;
+	}
+    	else{
+        	tol = 0.003;
+	}
+
+	double dArea1 = 0.0;
+
+	int nf = faces.size();
+
+    	for (int i = 0; i < nf; i++){
+
+	int nv_surf = faces[i].edges.size();
+		int nTriangles = faces[i].edges.size() - 2;
+		double h = MaxEdgeLength(faces[i],nodes);
+
+		if(h > tol){
+
+			FaceVector surf_fid;
+			NodeVector pnts_vor;
+			Node node;
+
+			for (int j = 0; j < nv_surf; j++){
+				pnts_vor.push_back(nodes[faces[i][j]]);
+			}
+
+	           	// insert elements and points
+            		int index = nv_surf;
+           		node = ((pnts_vor[0] + pnts_vor[1])/2.0)/(((pnts_vor[0] + pnts_vor[1])/2.0).Magnitude());
+            		pnts_vor.push_back(node);
+
+            		for (int j = 1; j < nv_surf-1; j++){
+
+				// insert elements
+				DataArray2D<int> surf_index;
+				surf_index.Allocate(4,3);
+
+				Face face1(3);
+				Face face2(3);
+				Face face3(3);
+				Face face4(3);
+
+				face1.SetNode(0,0), face1.SetNode(1,index), face1.SetNode(2,index+2);
+				face2.SetNode(0,index+2), face2.SetNode(1,index), face2.SetNode(2,index+1);
+				face3.SetNode(0,index+1), face3.SetNode(1,index), face3.SetNode(2,j);
+				face4.SetNode(0,index+2), face4.SetNode(1,index+1), face4.SetNode(2,j+1);
+
+				surf_fid.push_back(face1);
+				surf_fid.push_back(face2);
+				surf_fid.push_back(face3);
+				surf_fid.push_back(face4);
+
+                		index += 1;
+                		node = ((pnts_vor[j] + pnts_vor[j+1])/2.0)/(((pnts_vor[j] + pnts_vor[j+1])/2.0).Magnitude());
+                		pnts_vor.push_back(node);
+                		index += 1;
+                		node = ((pnts_vor[0] + pnts_vor[j+1])/2.0)/(((pnts_vor[0] + pnts_vor[j+1])/2.0).Magnitude());
+                		pnts_vor.push_back(node);		
+
+			}
+
+            		double dArea2 = CalculateFaceAreaQuadratureSplit(surf_fid,pnts_vor,nOrder);
+            		dArea1 += dArea2;
+
+		}
+		else{
+
+			double dArea2 = CalculateFaceAreaQuadrature(faces[i],nodes,nOrder);
+			dArea1 += dArea2;
+		}
+	}
+
+	return dArea1;
+
+
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+double MaxEdgeLength(
+	const Face & face,
+	const NodeVector & nodes
+) {
+	int nv_surf = face.edges.size();
+
+    	double h = 0.0;
+
+    	for (int i = 0; i < nv_surf-1; i++){
+		Node node_diff = nodes[face[i]]-nodes[face[i+1]];
+		h = fmax(h,node_diff.Magnitude());
+	}
+	Node node_diff = nodes[face[0]]-nodes[face[nv_surf-1]];
+	h = fmax(h,node_diff.Magnitude());
+	return h;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
 Real CalculateFaceAreaKarneysMethod(
 	const Face & face,
 	const NodeVector & nodes

src/GridElements.h

@@ -1051,5 +1051,40 @@ inline Node InterpolateQuadrilateralNode(
 
 ///////////////////////////////////////////////////////////////////////////////
 
+///	<summary>
+///		Find the maximum edge length of an element.
+///	</summary>
+
+double MaxEdgeLength(
+	const Face & face,
+	const NodeVector & nodes
+);
+
+///////////////////////////////////////////////////////////////////////////////
+
+///	<summary>
+///		Calculates the area of an element; to be used in adaptive element area calculation.
+///	</summary>
+
+Real CalculateFaceAreaQuadrature(
+	const Face & face,
+	const NodeVector & nodes,
+	int nOrder
+);
+
+///////////////////////////////////////////////////////////////////////////////
+
+///	<summary>
+///		Adaptive area calculation that refines an element based on the maximum edge length.
+///	</summary>
+
+Real CalculateFaceAreaQuadratureSplit(
+	const FaceVector & faces,
+	const NodeVector & nodes,
+	int & nOrder
+);
+
+///////////////////////////////////////////////////////////////////////////////
+
 #endif