diff --git a/src/GaussQuadrature.cpp b/src/GaussQuadrature.cpp index ab51d6e..4fcfd05 100644 --- a/src/GaussQuadrature.cpp +++ b/src/GaussQuadrature.cpp @@ -218,3 +218,63 @@ void GaussQuadrature::GetPoints( /////////////////////////////////////////////////////////////////////////////// +void GaussQuadrature::GetPoints( + int nCount, + DataArray2D & dG, + DataArray1D & 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; + + } + +} + diff --git a/src/GaussQuadrature.h b/src/GaussQuadrature.h index c743a12..59e567d 100644 --- a/src/GaussQuadrature.h +++ b/src/GaussQuadrature.h @@ -18,6 +18,7 @@ #define _GAUSSQUADRATURE_H_ #include "DataArray1D.h" +#include "DataArray2D.h" /////////////////////////////////////////////////////////////////////////////// @@ -48,6 +49,17 @@ class GaussQuadrature { DataArray1D & dG, DataArray1D & dW ); + + /// + /// Return the two-dimensional Gauss quadrature points and their corresponding + /// weights for the given number of points and reference element. + /// + + static void GetPoints( + int nCount, + DataArray2D & dG, + DataArray1D & dW + ); }; /////////////////////////////////////////////////////////////////////////////// diff --git a/src/GridElements.cpp b/src/GridElements.cpp index 60e8883..bfcbc78 100644 --- a/src/GridElements.cpp +++ b/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 dG1; + DataArray1D 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 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 diff --git a/src/GridElements.h b/src/GridElements.h index 3cd4fe5..1336c40 100644 --- a/src/GridElements.h +++ b/src/GridElements.h @@ -1051,5 +1051,40 @@ inline Node InterpolateQuadrilateralNode( /////////////////////////////////////////////////////////////////////////////// +/// +/// Find the maximum edge length of an element. +/// + +double MaxEdgeLength( + const Face & face, + const NodeVector & nodes +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Calculates the area of an element; to be used in adaptive element area calculation. +/// + +Real CalculateFaceAreaQuadrature( + const Face & face, + const NodeVector & nodes, + int nOrder +); + +/////////////////////////////////////////////////////////////////////////////// + +/// +/// Adaptive area calculation that refines an element based on the maximum edge length. +/// + +Real CalculateFaceAreaQuadratureSplit( + const FaceVector & faces, + const NodeVector & nodes, + int & nOrder +); + +/////////////////////////////////////////////////////////////////////////////// + #endif