You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 
 

65 lines
1.7 KiB

#include <igl/boundary_facets.h>
#include <igl/colon.h>
#include <igl/cotmatrix.h>
#include <igl/jet.h>
#include <igl/min_quad_with_fixed.h>
#include <igl/readOFF.h>
#include <igl/setdiff.h>
#include <igl/slice.h>
#include <igl/unique.h>
#include <igl/opengl/glfw/Viewer.h>
#include <Eigen/Sparse>
#include <iostream>
int main(int argc, char *argv[])
{
using namespace Eigen;
using namespace std;
MatrixXd V;
MatrixXi F;
igl::readOFF(TUTORIAL_SHARED_PATH "/camelhead.off",V,F);
// Find boundary edges
MatrixXi E;
igl::boundary_facets(F,E);
// Find boundary vertices
VectorXi b,IA,IC;
igl::unique(E,b,IA,IC);
// List of all vertex indices
VectorXi all,in;
igl::colon<int>(0,V.rows()-1,all);
// List of interior indices
igl::setdiff(all,b,in,IA);
// Construct and slice up Laplacian
SparseMatrix<double> L,L_in_in,L_in_b;
igl::cotmatrix(V,F,L);
igl::slice(L,in,in,L_in_in);
igl::slice(L,in,b,L_in_b);
// Dirichlet boundary conditions from z-coordinate
VectorXd Z = V.col(2);
VectorXd bc = Z(b);
// Solve PDE
SimplicialLLT<SparseMatrix<double > > solver(-L_in_in);
// slice into solution
Z(in) = solver.solve(L_in_b*bc);
// Alternative, short hand
igl::min_quad_with_fixed_data<double> mqwf;
// Linear term is 0
VectorXd B = VectorXd::Zero(V.rows(),1);
// Empty constraints
VectorXd Beq;
SparseMatrix<double> Aeq;
// Our cotmatrix is _negative_ definite, so flip sign
igl::min_quad_with_fixed_precompute((-L).eval(),b,Aeq,true,mqwf);
igl::min_quad_with_fixed_solve(mqwf,B,bc,Beq,Z);
// Plot the mesh with pseudocolors
igl::opengl::glfw::Viewer viewer;
viewer.data().set_mesh(V, F);
viewer.data().show_lines = false;
viewer.data().set_data(Z);
viewer.launch();
}