| 14 | |
| 15 | |
| 16 | int main() |
| 17 | { |
| 18 | bool dirichletBoundary = false; |
| 19 | directional::readOFF(TUTORIAL_DATA_PATH "/flatmesh.off",mesh); |
| 20 | |
| 21 | Eigen::SparseMatrix<double> d0 = directional::d0_matrix<double>(mesh, false); |
| 22 | Eigen::SparseMatrix<double> d1 = directional::d1_matrix<double>(mesh, false); |
| 23 | Eigen::SparseMatrix<double> h1, invh1, h2, invh2; |
| 24 | directional::hodge_star_1_matrix(mesh, h1, invh1); |
| 25 | directional::hodge_star_2_matrix<double>(mesh, h2, invh2); |
| 26 | |
| 27 | //demonstrating the exact sequences even with boundary conditions |
| 28 | Eigen::SparseMatrix<double> d1d0 = d1*d0; |
| 29 | double maxAbsValue = 0.0; |
| 30 | for (int k = 0; k < d1d0.outerSize(); ++k) |
| 31 | for (Eigen::SparseMatrix<double>::InnerIterator it(d1d0, k); it; ++it) |
| 32 | maxAbsValue = std::max(maxAbsValue, std::abs(it.value())); |
| 33 | |
| 34 | std::cout<<"exact sequence identity (should be exactly zero): "<<maxAbsValue<<std::endl; |
| 35 | |
| 36 | Eigen::MatrixXd harmBasis; |
| 37 | int bettiNumber = mesh.EV.rows() - (mesh.V.rows()-1) - mesh.F.rows(); |
| 38 | directional::cohomology_basis(d0, d1, h1, bettiNumber, harmBasis); |
| 39 | std::cout<<"Betti number: "<<bettiNumber<<std::endl; |
| 40 | |
| 41 | Eigen::RowVector3d COM = mesh.V.colwise().mean(); |
| 42 | Eigen::VectorXd z0GT(mesh.dcel.vertices.size()), curlGT(mesh.dcel.faces.size()); |
| 43 | for (int i=0;i<mesh.dcel.vertices.size();i++) |
| 44 | z0GT[i] = 3.0*cos((mesh.V(i,0)-COM(1))/4.0)*cos((mesh.V(i,2)-COM(2))/4.0); |
| 45 | |
| 46 | //vertexVec.array()/=vertexVec.mean(); |
| 47 | |
| 48 | for (int i=0;i<mesh.dcel.faces.size();i++){ |
| 49 | Eigen::RowVector3d midFacePointCOM = mesh.barycenters.row(i) - COM; |
| 50 | curlGT[i] = mesh.faceAreas[i]*(sin(midFacePointCOM(0)/2.0)+sin(midFacePointCOM(1)/2.0)+sin(midFacePointCOM(2)/2.0)); |
| 51 | } |
| 52 | curlGT.array()-=curlGT.mean(); //Due to Neumann boundary conditions (tangent coexact field), vector potential function adds up to zero |
| 53 | |
| 54 | //Generating an artificial composition and then reproducing it through decomposition. |
| 55 | Eigen::VectorXd harmGT = harmBasis.col(0); |
| 56 | Eigen::VectorXd z1ExactGT = d0*z0GT; |
| 57 | Eigen::VectorXd z1CoexactGT, curlFiltered; |
| 58 | directional::project_exact(d1, h1, curlGT, z1CoexactGT, curlFiltered, true); |
| 59 | |
| 60 | //creating balanced GT results for exact, coexact, and harmonic |
| 61 | z1ExactGT.array()/=sqrt(((z1ExactGT.transpose()*h1*z1ExactGT).coeff(0,0))); |
| 62 | z1CoexactGT.array()/=sqrt(((z1CoexactGT.transpose()*h1*z1CoexactGT).coeff(0,0))); |
| 63 | |
| 64 | z1 = harmGT + z1ExactGT + z1CoexactGT; |
| 65 | |
| 66 | Eigen::SparseMatrix<double> I(mesh.F.rows(), mesh.F.rows()); |
| 67 | I.setIdentity(); |
| 68 | directional::hodge_decomposition<double>(d0, d1, h1, h2, z1, 1, z1Exact, z1Coexact, z1Harmonic, z0, z2); |
| 69 | std::cout<<"Exact reproduction: "<<(z1Exact - z1ExactGT).cwiseAbs().maxCoeff()<<std::endl; |
| 70 | std::cout<<"Coexact reproduction: "<<(z1Coexact - z1CoexactGT).cwiseAbs().maxCoeff()<<std::endl; |
| 71 | std::cout<<"Harmonic reproduction: "<<(z1Harmonic - harmGT).cwiseAbs().maxCoeff()<<std::endl; |
| 72 | |
| 73 | //triangle mesh setup |
nothing calls this directly
no test coverage detected