Corrected calculus

ChangeLog.md

@@ -1,3 +1,10 @@
+# DGtal 1.5beta
+## New features / critical changes
+## Changes
+- *DEC*
+  - Add normal corrected FEM calculus for Poisson problem and Interpolated Calculus variant of DEC (Colin Weill--Duflos, PULL_REQUEST)
+## Bug fixes
+
 # DGtal 1.4
 
 ## New features / critical changes
@@ -8,8 +15,8 @@
     conan.io, especially on Windows, new `ENABLE_CONAN` cmake
     option to activate this. (David Coeurjolly,
     [#1689](https://github.com/DGtal-team/DGtal/pull/1689))
-  - Faster build using CPM for dependency download and ccache with the cmake `USE_CCACHE=YES`option 
-    (ccache must be installed). (David Coeurjolly, [#1696](https://github.com/DGtal-team/DGtal/pull/1696)) 
+  - Faster build using CPM for dependency download and ccache with the cmake `USE_CCACHE=YES`option
+    (ccache must be installed). (David Coeurjolly, [#1696](https://github.com/DGtal-team/DGtal/pull/1696))
   - Better documentation style using doxygen-awesome.css. (David Coeurjolly,
     [#1697](https://github.com/DGtal-team/DGtal/pull/1697))
 
@@ -25,8 +32,8 @@
     [#1721](https://github.com/DGtal-team/DGtal/pull/1721))
   - Add CMake option DGTAL_WRAP_PYTHON (Pablo Hernandez-Cerdan,
     [#1700](https://github.com/DGtal-team/DGtal/pull/1700))
-  - Upgrade of the conan scripts (for windows build) to conan 2, removing the ENABLE_CONAN option 
-    (documentation update instead) (David Coeurjolly, 
+  - Upgrade of the conan scripts (for windows build) to conan 2, removing the ENABLE_CONAN option
+    (documentation update instead) (David Coeurjolly,
     [#1729](https://github.com/DGtal-team/DGtal/pull/1729))
 
 - *IO*
@@ -76,8 +83,8 @@
     Hull files (David Coeurjolly, [#1716](https://github.com/DGtal-team/DGtal/pull/1716))
   - Activate and fix CTest tests on windows system. (Bertrand Kerautret,
     [#1706](https://github.com/DGtal-team/DGtal/pull/1706))
-  - For now, removing Cairo deps install on windows (6hours long build 
-    with conan in the windows debug mode). (David Coeurjolly, 
+  - For now, removing Cairo deps install on windows (6hours long build
+    with conan in the windows debug mode). (David Coeurjolly,
     [#1705](https://github.com/DGtal-team/DGtal/pull/1705))
   - Fix conan file upload issue and log message. (Bertrand Kerautret,
     [#1704](https://github.com/DGtal-team/DGtal/pull/1704))
@@ -115,16 +122,16 @@
     [#1717](https://github.com/DGtal-team/DGtal/pull/1717))
   - Fix const attribute that shouldn't be in FreemanChain (Colin Weill--Duflos,
 	  [#1723](https://github.com/DGtal-team/DGtal/pull/1723))
-  - Fix seg fault due to recent compilers in FrechetShortcut (Bertrand Kerautret, 
+  - Fix seg fault due to recent compilers in FrechetShortcut (Bertrand Kerautret,
      Isabelle Sivignon [#1726](https://github.com/DGtal-team/DGtal/pull/1726))
-  - Fix FrechetShortcut to enable the parameter error to be equal to 0 and add new 
+  - Fix FrechetShortcut to enable the parameter error to be equal to 0 and add new
     tests in testFrechetShortcut (Isabelle Sivignon, [#1726](https://github.com/DGtal-team/DGtal/pull/1726))
 
 
 - *IO*
   - Fix of the `getHSV` method in the `Color` class. (David Coeurjolly,
     [#1674](https://github.com/DGtal-team/DGtal/pull/1674))
-  - Fix of `SurfaceMeshWriter::writeIsoLinesOBJ` 
+  - Fix of `SurfaceMeshWriter::writeIsoLinesOBJ`
     (Jacques-Olivier Lachaud, [#1701](https://github.com/DGtal-team/DGtal/pull/1701))
   - Fix of the `PointListReader::getPolygonsFromInputStream` (Xun Gong,
     [#1708](https://github.com/DGtal-team/DGtal/pull/1708))
@@ -538,7 +545,7 @@
     (Jacques-Olivier Lachaud,[#1470](https://github.com/DGtal-team/DGtal/pull/1470))
 
 - *io*
-  - The GenericWriter can now export in 3D ITK format (nii, mha,  mhd,  tiff).  
+  - The GenericWriter can now export in 3D ITK format (nii, mha,  mhd,  tiff).
     (Bertrand Kerautret [#1485](https://github.com/DGtal-team/DGtal/pull/1485))
   - New Viridis ColorGradientPreset and clean of  useless template specializations in
     the GenericWriter for color image. (Bertrand Kerautret
@@ -714,7 +721,7 @@
   - Fix apple clang  compilation issue with a workaround to the
     ConstIteratorAdapter class that does not satisfy the _is_forward concept of the STL:
     using boost::first_max_element instead std::max_element.
-    (Bertrand Kerautret, [#1437](https://github.com/DGtal-team/DGtal/pull/1437))  
+    (Bertrand Kerautret, [#1437](https://github.com/DGtal-team/DGtal/pull/1437))
   - Abort compilation at configure time when the compiler is gcc 10.1 due to compiler bug.
     Fix issue #1501.
     (Pablo Hernandez-Cerdan, [#1506](https://github.com/DGtal-team/DGtal/pull/1506))

examples/polyscope-examples/CMakeLists.txt

@@ -39,11 +39,13 @@ set(POLYSCOPE_SRCS
   exampleBunnyHead
   tangency-explorer
   tangency-reconstruction
+  dgtalFEM-poisson
+  dgtalCC-poisson
 )
 
 if(WITH_LIBIGL)
- set(POLYSCOPE_SRCS  
-      ${POLYSCOPE_SRCS} 
+ set(POLYSCOPE_SRCS
+      ${POLYSCOPE_SRCS}
       windingNumbersShape
  )
 endif()

examples/polyscope-examples/dgtalCC-poisson.cpp

@@ -0,0 +1,137 @@
+/**
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Lesser General Public License as
+ *  published by the Free Software Foundation, either version 3 of the
+ *  License, or  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+/**
+ * @file
+ * @author Colin Weill--Duflos (\c [email protected] )
+ * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
+ *
+ * @date 2024/06/05
+ *
+ * This file is part of the DGtal library.
+ */
+#include <iostream>
+#include <DGtal/base/Common.h>
+#include <DGtal/helpers/StdDefs.h>
+#include <DGtal/helpers/Shortcuts.h>
+#include <DGtal/helpers/ShortcutsGeometry.h>
+#include <DGtal/shapes/SurfaceMesh.h>
+#include <DGtal/geometry/surfaces/DigitalSurfaceRegularization.h>
+#include <DGtal/dec/InterpolatedCorrectedCalculus.h>
+#include <DGtal/math/linalg/DirichletConditions.h>
+
+#include <polyscope/polyscope.h>
+#include <polyscope/surface_mesh.h>
+#include <polyscope/point_cloud.h>
+
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+using namespace DGtal;
+using namespace Z3i;
+
+// Using standard 3D digital space.
+typedef Shortcuts<Z3i::KSpace>         SH3;
+typedef ShortcutsGeometry<Z3i::KSpace> SHG3;
+// The following typedefs are useful
+typedef SurfaceMesh< RealPoint, RealVector >  SurfMesh;
+typedef std::size_t             Index;
+//Polyscope global
+polyscope::SurfaceMesh *psMesh;
+SurfMesh surfmesh;
+float scale = 0.1;
+
+void computeLaplace()
+{
+  //! [CC-init]
+  typedef InterpolatedCorrectedCalculus<EigenLinearAlgebraBackend, SH3::RealPoint,SH3::RealVector> CC;
+  typedef DirichletConditions< EigenLinearAlgebraBackend >  DC;
+  CC calculus(surfmesh);
+  CC::LinearOperator L = calculus.L0();
+  CC::DenseVector g = CC::DenseVector::Zero(surfmesh.nbVertices());
+  DC::IntegerVector b = DC::IntegerVector::Zero( g.rows() );
+
+  //We set values on the boundary
+  auto boundaryEdges = surfmesh.computeManifoldBoundaryEdges();
+  std::cout<< "Number of boundary edges= "<<boundaryEdges.size()<<std::endl;
+
+  auto pihVertex=[&](const SurfMesh::Vertex &v){return  cos(scale*(surfmesh.position(v)[0]))*(scale*surfmesh.position(v)[1]);};
+
+  for(auto &e: boundaryEdges)
+  {
+    auto adjVertices = surfmesh.edgeVertices(e);
+    g(adjVertices.first)  = pihVertex(adjVertices.first);
+    g(adjVertices.second) = pihVertex(adjVertices.second);
+    b(adjVertices.first)  = 1;
+    b(adjVertices.second) = 1;
+  }
+
+  // Solve Δu=0 with g as boundary conditions
+  CC::LinearAlgebraBackend::SolverSimplicialLDLT solver;
+  CC::LinearOperator L_dirichlet = DC::dirichletOperator( L, b );
+  solver.compute( L_dirichlet );
+  ASSERT(solver.info()==Eigen::Success);
+  CC::DenseVector g_dirichlet = DC::dirichletVector( L, g, b, g );
+  CC::DenseVector x_dirichlet = solver.solve( g_dirichlet );
+  CC::DenseVector u = DC::dirichletSolution( x_dirichlet, b, g );
+  //! [CC-init]
+
+  psMesh->addVertexScalarQuantity("g", g);
+  psMesh->addVertexScalarQuantity("u", u);
+}
+
+void myCallback()
+{
+  ImGui::SliderFloat("Phi scale", &scale, 0., 10.);
+  if(ImGui::Button("Compute Laplace problem"))
+    computeLaplace();
+}
+
+int main()
+{
+  auto params = SH3::defaultParameters() | SHG3::defaultParameters() |  SHG3::parametersGeometryEstimation();
+
+  auto h=.7   ; //gridstep
+  params( "polynomial", "0.1*y*y -0.1*x*x - 2.0*z" )( "gridstep", h );
+  auto implicit_shape  = SH3::makeImplicitShape3D  ( params );
+  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
+  auto K               = SH3::getKSpace( params );
+  auto binary_image    = SH3::makeBinaryImage( digitized_shape, params );
+  auto surface         = SH3::makeDigitalSurface( binary_image, K, params );
+  auto primalSurface   = SH3::makePrimalSurfaceMesh(surface);
+
+  //Need to convert the faces
+  std::vector<std::vector<SH3::SurfaceMesh::Vertex>> faces;
+  std::vector<RealPoint> positions = primalSurface->positions();
+  for(auto face= 0 ; face < primalSurface->nbFaces(); ++face)
+    faces.push_back(primalSurface->incidentVertices( face ));
+
+  surfmesh = SurfMesh(positions.begin(),
+                      positions.end(),
+                      faces.begin(),
+                      faces.end());
+  surfmesh.computeFaceNormalsFromPositions();
+  surfmesh.computeVertexNormalsFromFaceNormals();
+  psMesh = polyscope::registerSurfaceMesh("digital surface", positions, faces);
+
+  // Initialize polyscope
+  polyscope::init();
+
+  // Set the callback function
+  polyscope::state::userCallback = myCallback;
+  polyscope::show();
+  return EXIT_SUCCESS;
+}

examples/polyscope-examples/dgtalFEM-poisson.cpp

@@ -0,0 +1,136 @@
+/**
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Lesser General Public License as
+ *  published by the Free Software Foundation, either version 3 of the
+ *  License, or  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+/**
+ * @file
+ * @author Colin Weill--Duflos (\c [email protected] )
+ * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
+ *
+ * @date 2024/06/05
+ *
+ * This file is part of the DGtal library.
+ */
+#include <iostream>
+#include <DGtal/base/Common.h>
+#include <DGtal/helpers/StdDefs.h>
+#include <DGtal/helpers/Shortcuts.h>
+#include <DGtal/helpers/ShortcutsGeometry.h>
+#include <DGtal/shapes/SurfaceMesh.h>
+#include <DGtal/geometry/surfaces/DigitalSurfaceRegularization.h>
+#include <DGtal/dec/NormalCorrectedFEM.h>
+#include <DGtal/math/linalg/DirichletConditions.h>
+
+#include <polyscope/polyscope.h>
+#include <polyscope/surface_mesh.h>
+#include <polyscope/point_cloud.h>
+
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+using namespace DGtal;
+using namespace Z3i;
+
+// Using standard 3D digital space.
+typedef Shortcuts<Z3i::KSpace>         SH3;
+typedef ShortcutsGeometry<Z3i::KSpace> SHG3;
+// The following typedefs are useful
+typedef SurfaceMesh< RealPoint, RealVector >  SurfMesh;
+typedef std::size_t             Index;
+//Polyscope global
+polyscope::SurfaceMesh *psMesh;
+SurfMesh surfmesh;
+float scale = 0.1;
+
+void computeLaplace()
+{
+  //! [FEM-init]
+  typedef NormalCorrectedFEM<EigenLinearAlgebraBackend, SH3::RealPoint,SH3::RealVector> ncFEM;
+  typedef DirichletConditions< EigenLinearAlgebraBackend >  DC;
+  ncFEM calculus(surfmesh);
+  ncFEM::LinearOperator L = calculus.L0();
+  ncFEM::DenseVector g = ncFEM::DenseVector::Zero(surfmesh.nbVertices());
+  DC::IntegerVector b = DC::IntegerVector::Zero( g.rows() );
+
+  //We set values on the boundary
+  auto boundaryEdges = surfmesh.computeManifoldBoundaryEdges();
+  std::cout<< "Number of boundary edges= "<<boundaryEdges.size()<<std::endl;
+
+  auto pihVertex=[&](const SurfMesh::Vertex &v){return  cos(scale*(surfmesh.position(v)[0]))*(scale*surfmesh.position(v)[1]);};
+
+  for(auto &e: boundaryEdges)
+  {
+    auto adjVertices = surfmesh.edgeVertices(e);
+    g(adjVertices.first)  = pihVertex(adjVertices.first);
+    g(adjVertices.second) = pihVertex(adjVertices.second);
+    b(adjVertices.first)  = 1;
+    b(adjVertices.second) = 1;
+  }
+
+  // Solve Δu=0 with g as boundary conditions
+  ncFEM::LinearAlgebraBackend::SolverSimplicialLDLT solver;
+  ncFEM::LinearOperator L_dirichlet = DC::dirichletOperator( L, b );
+  solver.compute( L_dirichlet );
+  ASSERT(solver.info()==Eigen::Success);
+  ncFEM::DenseVector g_dirichlet = DC::dirichletVector( L, g, b, g );
+  ncFEM::DenseVector x_dirichlet = solver.solve( g_dirichlet );
+  ncFEM::DenseVector u = DC::dirichletSolution( x_dirichlet, b, g );
+  //! [FEM-init]
+
+  psMesh->addVertexScalarQuantity("g", g);
+  psMesh->addVertexScalarQuantity("u", u);
+}
+
+void myCallback()
+{
+  ImGui::SliderFloat("Phi scale", &scale, 0., 10.);
+  if(ImGui::Button("Compute Laplace problem"))
+    computeLaplace();
+}
+
+int main()
+{
+  auto params = SH3::defaultParameters() | SHG3::defaultParameters() |  SHG3::parametersGeometryEstimation();
+
+  auto h=.7   ; //gridstep
+  params( "polynomial", "0.1*y*y -0.1*x*x - 2.0*z" )( "gridstep", h );
+  auto implicit_shape  = SH3::makeImplicitShape3D  ( params );
+  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
+  auto K               = SH3::getKSpace( params );
+  auto binary_image    = SH3::makeBinaryImage( digitized_shape, params );
+  auto surface         = SH3::makeDigitalSurface( binary_image, K, params );
+  auto primalSurface   = SH3::makePrimalSurfaceMesh(surface);
+
+  //Need to convert the faces
+  std::vector<std::vector<SH3::SurfaceMesh::Vertex>> faces;
+  std::vector<RealPoint> positions = primalSurface->positions();
+  for(auto face= 0 ; face < primalSurface->nbFaces(); ++face)
+    faces.push_back(primalSurface->incidentVertices( face ));
+
+  surfmesh = SurfMesh(positions.begin(),
+                      positions.end(),
+                      faces.begin(),
+                      faces.end());
+  surfmesh.computeFaceNormalsFromPositions();
+  psMesh = polyscope::registerSurfaceMesh("digital surface", positions, faces);
+
+  // Initialize polyscope
+  polyscope::init();
+
+  // Set the callback function
+  polyscope::state::userCallback = myCallback;
+  polyscope::show();
+  return EXIT_SUCCESS;
+}