diff --git a/geometry-processing/6-GeomProcessing/c++/CMakeLists.txt b/geometry-processing/6-GeomProcessing/c++/CMakeLists.txt index 604a626..61920b3 100644 --- a/geometry-processing/6-GeomProcessing/c++/CMakeLists.txt +++ b/geometry-processing/6-GeomProcessing/c++/CMakeLists.txt @@ -1,3 +1,6 @@ add_executable(geomproc geomproc.cpp) target_link_libraries(geomproc polyscope geometry-central) + +add_executable(dirichletProblems dirichletProblems.cpp) +target_link_libraries(dirichletProblems polyscope geometry-central) \ No newline at end of file diff --git a/geometry-processing/6-GeomProcessing/c++/dirichletProblems.cpp b/geometry-processing/6-GeomProcessing/c++/dirichletProblems.cpp new file mode 100644 index 0000000..670fbec --- /dev/null +++ b/geometry-processing/6-GeomProcessing/c++/dirichletProblems.cpp @@ -0,0 +1,132 @@ +#include +#include +#include + +#include +#include + +using namespace geometrycentral; +using namespace geometrycentral::surface; + +// == Geometry-central data +std::unique_ptr mesh; +std::unique_ptr geometry; + +/** Gradient of fvals, a scalar quantity on vertices. Outputs a vector quantity on one face */ +Vector3 face_gradient_of(Face f, std::vector fvals) { + Vector3 normal = geometry->faceNormal(f); + double area = geometry->faceArea(f); + + std::vector v_indices; + for(Vertex v: f.adjacentVertices()) { + v_indices.push_back(v); + } + + Vector3 loc_face_gradient = {0, 0, 0}; + for (int i=0; i < 3; i++) { + Vertex v0 = v_indices.at(i); + Vertex v1 = v_indices.at((i+1)%3); + Vertex v2 = v_indices.at((i+2)%3); + + loc_face_gradient += ( + fvals[v0.getIndex()]/(2*area)* + cross(normal, geometry->vertexPositions[v2]-geometry->vertexPositions[v1]) + ); + } + return loc_face_gradient; +} + +/** Compute the divergence of U (vector on faces), outputs as result a scalar vertices */ +std::vector divergence(std::vector U) { + std::vector out; + + for (Vertex v : mesh->vertices()) { + out.push_back(0); + } + + for (Face f : mesh->faces()) { + Vector3 normal = geometry->faceNormal(f); + std::vector v_indices; + for(Vertex v: f.adjacentVertices()) { + v_indices.push_back(v); + } + + for (int i=0; i < 3; i++) { + Vertex v0 = v_indices.at(i); + Vertex v1 = v_indices.at((i+1)%3); + Vertex v2 = v_indices.at((i+2)%3); + + out[v0.getIndex()] -= dot( + U[f.getIndex()], + cross(normal, geometry->vertexPositions[v2]-geometry->vertexPositions[v1]) + ); + } + } + return out; +} + +int main(int argc, char **argv) { + // Initialize polyscope + polyscope::init(); + + // Load mesh + std::tie(mesh, geometry) = readManifoldSurfaceMesh(argv[1]); + + // Register the mesh with polyscope + auto sfm = polyscope::registerSurfaceMesh("Input obj", + geometry->inputVertexPositions, + mesh->getFaceVertexList(), + polyscopePermutations(*mesh)); + geometry->requireVertexPositions(); + geometry->requireFaceNormals(); + geometry->requireFaceAreas(); + + std::vector sin_coords; + {//* First question : showing basic quantity + std::vector snd_coordinate; + for (auto i=0; i < geometry->inputVertexPositions.size(); i++) { + auto pos = geometry->inputVertexPositions[i]; + snd_coordinate.push_back((double)pos[1]); + } + sfm->addVertexScalarQuantity("Second coordinate", snd_coordinate); + + float lambda = 18.88; + for (auto i=0; i < geometry->inputVertexPositions.size(); i++) { + auto pos = geometry->inputVertexPositions[i]; + double posX = cos(PI / 4) * pos.x - sin(PI / 4) * pos.y; + double posY = sin(PI / 4) * pos.x + cos(PI / 4) * pos.y; + sin_coords.push_back(sin(lambda * (posX + 0.5)) + sin(lambda * posY)); + } + sfm->addVertexScalarQuantity("Sinusoid data from coordinates", sin_coords); + } + + {//* Second question : compute the area and the normal vector for each triangle + // After doing geometry->requireFaceNormals(), + // we can access normals with geometry->faceNormal(Face f) + + sfm->addFaceVectorQuantity("Faces normals", geometry->faceNormals); + + // After doing geometry->requireFaceAreas(), + // we can access normals with geometry->faceArea(Face f) + sfm->addFaceScalarQuantity("Faces area", geometry->faceAreas); + } + + std::vector face_gradients; + {//* Third question : compute the per face gradient of f + for (Face f: mesh->faces()) { + face_gradients.push_back(face_gradient_of(f, sin_coords)); + } + sfm->addFaceVectorQuantity("Faces gradients", face_gradients); + } + + {//* Fourth question : Compute and display the divergence of grad(f) + sfm->addVertexScalarQuantity("Faces div of gradients", divergence(face_gradients)); + } + + //* Laplace-Beltrami + + // Give control to the polyscope gui + polyscope::show(); + + return EXIT_SUCCESS; +} diff --git a/geometry-processing/6-GeomProcessing/c++/geomproc.cpp b/geometry-processing/6-GeomProcessing/c++/geomproc.cpp index 0c8ccb3..42a3bb7 100644 --- a/geometry-processing/6-GeomProcessing/c++/geomproc.cpp +++ b/geometry-processing/6-GeomProcessing/c++/geomproc.cpp @@ -22,13 +22,28 @@ int main(int argc, char **argv) { std::tie(mesh, geometry) = readManifoldSurfaceMesh(argv[1]); // Register the mesh with polyscope - polyscope::registerSurfaceMesh("Input obj", + auto pcl = polyscope::registerSurfaceMesh("Input obj", geometry->inputVertexPositions, mesh->getFaceVertexList(), polyscopePermutations(*mesh)); + geometry->requireVertexPositions(); + + std::vector nbAdjFaces; + for (auto i=0; i < geometry->inputVertexPositions.size(); i++) { + nbAdjFaces.push_back(0); + } - - + for(Face f: mesh->faces()) { + std::cout << " Je suis sur la face " << f.getIndex() << " : " << std::endl; + for(Vertex v: f.adjacentVertices()) { + auto coord = geometry->vertexPositions[v]; + std::cout << " vertex: "<addVertexScalarQuantity("nb of adjacent faces", nbAdjFaces); // Give control to the polyscope gui polyscope::show();