boundaryedges and CR mass matrix

Former-commit-id: c3c051a610b770dd9314bd1c7327f29812d7b574

include/igl/crouzeix_raviart_massmatrix.cpp

@@ -0,0 +1,69 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "crouzeix_raviart_massmatrix.h"
+#include "unique_simplices.h"
+#include "all_edges.h"
+
+#include "is_edge_manifold.h"
+#include "doublearea.h"
+
+#include <cassert>
+#include <vector>
+
+template <typename MT, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedEMAP>
+void igl::crouzeix_raviart_massmatrix(
+    const Eigen::PlainObjectBase<DerivedV> & V, 
+    const Eigen::PlainObjectBase<DerivedF> & F, 
+    Eigen::SparseMatrix<MT> & M,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
+{
+  // All occurances of directed edges
+  Eigen::MatrixXi allE;
+  all_edges(F,allE);
+  Eigen::VectorXi _1;
+  unique_simplices(allE,E,_1,EMAP);
+  return crouzeix_raviart_massmatrix(V,F,E,EMAP,M);
+}
+
+template <typename MT, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedEMAP>
+void igl::crouzeix_raviart_massmatrix(
+    const Eigen::PlainObjectBase<DerivedV> & V, 
+    const Eigen::PlainObjectBase<DerivedF> & F, 
+    const Eigen::PlainObjectBase<DerivedE> & E,
+    const Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+    Eigen::SparseMatrix<MT> & M)
+{
+  using namespace igl;
+  using namespace Eigen;
+  using namespace std;
+  assert(F.cols() == 3);
+  // Mesh should be edge-manifold
+  assert(is_edge_manifold(V,F));
+  // number of elements (triangles)
+  int m = F.rows();
+  // Get triangle areas
+  VectorXd TA;
+  doublearea(V,F,TA);
+  vector<Triplet<MT> > MIJV(3*m);
+  for(int f = 0;f<m;f++)
+  {
+    for(int c = 0;c<3;c++)
+    {
+      MIJV[f+m*c] = Triplet<MT>(EMAP(f+m*c),EMAP(f+m*c),TA(f)*0.5);
+    }
+  }
+  M.resize(E.rows(),E.rows());
+  M.setFromTriplets(MIJV.begin(),MIJV.end());
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit instanciation
+template void igl::crouzeix_raviart_massmatrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::crouzeix_raviart_massmatrix<float, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<float, 0, int>&);
+#endif

include/igl/crouzeix_raviart_massmatrix.h

@@ -0,0 +1,60 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef CROUZEIX_RAVIART_MASSMATRIX_H
+#define CROUZEIX_RAVIART_MASSMATRIX_H
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace igl
+{
+  // CROUZEIX_RAVIART_MASSMATRIX Compute the Crouzeix-Raviart mass matrix where
+  // M(e,e) is just the sum of the areas of the triangles on either side of an
+  // edge e.
+  //
+  // See for example "Discrete Quadratic Curvature Energies" [Wardetzky, Bergou,
+  // Harmon, Zorin, Grinspun 2007]
+  //
+  // Templates:
+  //   MT  type of eigen sparse matrix for M (e.g. double for
+  //     SparseMatrix<double>)
+  //   DerivedV  derived type of eigen matrix for V (e.g. derived from
+  //     MatrixXd)
+  //   DerivedF  derived type of eigen matrix for F (e.g. derived from
+  //     MatrixXi)
+  //   DerivedE  derived type of eigen matrix for E (e.g. derived from
+  //     MatrixXi)
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   F  #F by 3 list of triangle indices
+  // Outputs:
+  //   M  #E by #E edge-based diagonal mass matrix
+  //   E  #E by 2 list of edges
+  //   EMAP  #F*3 list of indices mapping allE to E
+  //
+  //
+  template <typename MT, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedEMAP>
+  void crouzeix_raviart_massmatrix(
+      const Eigen::PlainObjectBase<DerivedV> & V, 
+      const Eigen::PlainObjectBase<DerivedF> & F, 
+      Eigen::SparseMatrix<MT> & M,
+      Eigen::PlainObjectBase<DerivedE> & E,
+      Eigen::PlainObjectBase<DerivedEMAP> & EMAP);
+  // wrapper if E and EMAP are already computed (better match!)
+  template <typename MT, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedEMAP>
+  void crouzeix_raviart_massmatrix(
+      const Eigen::PlainObjectBase<DerivedV> & V, 
+      const Eigen::PlainObjectBase<DerivedF> & F, 
+      const Eigen::PlainObjectBase<DerivedE> & E,
+      const Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+      Eigen::SparseMatrix<MT> & M);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "crouzeix_raviart_massmatrix.cpp"
+#endif
+  
+#endif

include/igl/is_boundary_edge.cpp

@@ -0,0 +1,124 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "is_boundary_edge.h"
+#include "unique.h"
+#include "sort.h"
+
+template <
+  typename DerivedF,
+  typename DerivedE,
+  typename DerivedB>
+void igl::is_boundary_edge(
+  const Eigen::PlainObjectBase<DerivedE> & E,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedB> & B)
+{
+  using namespace igl;
+  using namespace Eigen;
+  using namespace std;
+  // Should be triangles
+  assert(F.cols() == 3);
+  // Should be edges
+  assert(E.cols() == 2);
+  // number of faces
+  const int m = F.rows();
+  // Collect all directed edges after E
+  MatrixXi EallE(E.rows()+3*m,2);
+  EallE.block(0,0,E.rows(),E.cols()) = E;
+  for(int e = 0;e<3;e++)
+  {
+    for(int f = 0;f<m;f++)
+    {
+      for(int c = 0;c<2;c++)
+      {
+        // 12 20 01
+        EallE(E.rows()+m*e+f,c) = F(f,(c+1+e)%3);
+      }
+    }
+  }
+  // sort directed edges into undirected edges
+  MatrixXi sorted_EallE,_;
+  sort(EallE,2,true,sorted_EallE,_);
+  // Determine unique undirected edges E and map to directed edges EMAP
+  MatrixXi uE;
+  VectorXi EMAP;
+  unique_rows(sorted_EallE,uE,_,EMAP);
+  // Counts of occurances
+  VectorXi N = VectorXi::Zero(uE.rows());
+  for(int e = 0;e<EMAP.rows();e++)
+  {
+    N(EMAP(e))++;
+  }
+  B.resize(E.rows());
+  // Look of occurances of 2: one for original and another for boundary
+  for(int e = 0;e<E.rows();e++)
+  {
+    B(e) = (N(EMAP(e)) == 2);
+  }
+}
+
+
+template <
+  typename DerivedF,
+  typename DerivedE,
+  typename DerivedB,
+  typename DerivedEMAP>
+void igl::is_boundary_edge(
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedB> & B,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
+{
+  using namespace igl;
+  using namespace Eigen;
+  using namespace std;
+  // Should be triangles
+  assert(F.cols() == 3);
+  // number of faces
+  const int m = F.rows();
+  // Collect all directed edges after E
+  MatrixXi allE(3*m,2);
+  for(int e = 0;e<3;e++)
+  {
+    for(int f = 0;f<m;f++)
+    {
+      for(int c = 0;c<2;c++)
+      {
+        // 12 20 01
+        allE(m*e+f,c) = F(f,(c+1+e)%3);
+      }
+    }
+  }
+  // sort directed edges into undirected edges
+  MatrixXi sorted_allE,_;
+  sort(allE,2,true,sorted_allE,_);
+  // Determine unique undirected edges E and map to directed edges EMAP
+  unique_rows(sorted_allE,E,_,EMAP);
+  // Counts of occurances
+  VectorXi N = VectorXi::Zero(E.rows());
+  for(int e = 0;e<EMAP.rows();e++)
+  {
+    N(EMAP(e))++;
+  }
+  B.resize(E.rows());
+  // Look of occurances of 1
+  for(int e = 0;e<E.rows();e++)
+  {
+    B(e) = (N(EMAP(e)) == 1);
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit instanciation:
+template void igl::is_boundary_edge<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::is_boundary_edge<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
+template void igl::is_boundary_edge<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
+template void igl::is_boundary_edge<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::is_boundary_edge<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::is_boundary_edge<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/is_boundary_edge.h

@@ -0,0 +1,51 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IS_BOUNDARY_EDGE_H
+#define IS_BOUNDARY_EDGE_H
+#include <Eigen/Dense>
+
+namespace igl
+{
+  //  IS_BOUNDARY_EDGE Determine for each edge E if it is a "boundary edge" in F.
+  //  Boundary edges are undirected edges which occur only once.
+  // 
+  //  Inputs:
+  //    E  #E by 2 list of edges
+  //    F  #F by 3 list of triangles
+  //  Outputs:
+  //    B  #E list bools. true iff unoriented edge occurs exactly once in F
+  //      (non-manifold and non-existant edges will be false)
+  // 
+  template <
+    typename DerivedF,
+    typename DerivedE,
+    typename DerivedB>
+  void is_boundary_edge(
+    const Eigen::PlainObjectBase<DerivedE> & E,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedB> & B);
+  // Wrapper where Edges should also be computed from F
+  //   E  #E by 2 list of edges
+  //   EMAP  #F*3 list of indices mapping allE to E
+  template <
+    typename DerivedF,
+    typename DerivedE,
+    typename DerivedB,
+    typename DerivedEMAP>
+  void is_boundary_edge(
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedB> & B,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "is_boundary_edge.cpp"
+#endif
+
+#endif