add code

CMakeLists.txt

@@ -0,0 +1,24 @@
+cmake_minimum_required(VERSION 3.16)
+project(ImplictSim LANGUAGES CXX)
+
+
+
+file(GLOB SRC_FILES 
+            src/BgMesh/*.cpp
+            src/Geometry/Implict/*.cpp
+            src/Quadrature/*.cpp
+            )
+
+add_library(${PROJECT_NAME}_dev ${SRC_FILES})
+
+target_include_directories(${PROJECT_NAME}_dev PUBLIC /usr/local/include)
+target_include_directories(${PROJECT_NAME}_dev PUBLIC include)
+target_include_directories(${PROJECT_NAME}_dev PUBLIC include/BgMesh)
+target_include_directories(${PROJECT_NAME}_dev PUBLIC include/Geometry)
+target_include_directories(${PROJECT_NAME}_dev PUBLIC include/Utils)
+
+add_executable(${PROJECT_NAME}_bin src/main.cpp)
+target_link_libraries(${PROJECT_NAME}_bin PUBLIC ${PROJECT_NAME}_dev)
+
+
+

include/BgMesh/BgMesh.hpp

@@ -0,0 +1,37 @@
+#pragma once
+
+#include <Eigen/Eigen>
+#include <vector>
+
+namespace ImplictSIM {
+
+template <size_t dim>
+class BgMeshData {
+public:
+    Eigen::Vector<double, dim> leftDownPnt;
+    Eigen::Vector<double, dim> centerPnt;
+    Eigen::Vector<double, dim> rightUpPnt;
+
+public:
+    BgMeshData(Eigen::Vector<double, dim> leftDownPnt_, Eigen::Vector<double, dim> rightUpPnt_)
+    {
+        leftDownPnt = leftDownPnt_;
+        rightUpPnt = rightUpPnt_;
+        centerPnt = (leftDownPnt + rightUpPnt) / 2.0;
+    }
+};
+
+template <size_t dim>
+class BgMesh {
+public:
+    int leafNum;
+    Eigen::Vector<double, dim> leftDownPnt;
+    Eigen::Vector<double, dim> rightUpPnt;
+    std::vector<BgMeshData<dim>> meshData;
+
+public:
+    BgMesh(Eigen::Vector<double, dim> leftDownPnt_, Eigen::Vector<double, dim> rightUpPnt_)
+    {
+    }
+};
+};

include/Geometry/Implict/ImplictModel.hpp

@@ -0,0 +1,47 @@
+#pragma once
+
+#include <Eigen/Eigen>
+
+namespace ImplictSIM {
+
+template <size_t dim>
+class ImplictModel {
+public:
+    typedef double (*disCalculationPtr)(const Eigen::Vector<double, dim>& point);
+    typedef double (*gradientDirXPtr)(const Eigen::Vector<double, dim>& point);
+    typedef double (*gradientDirYPtr)(const Eigen::Vector<double, dim>& point);
+    typedef double (*gradientDirZPtr)(const Eigen::Vector<double, dim>& point);
+    ImplictModel(disCalculationPtr disFunc_, gradientDirXPtr gradOnX_,
+        gradientDirYPtr gradOnY_, gradientDirZPtr gradOnZ_)
+        : disFunc(disFunc_)
+        , gradOnX(gradOnX_)
+        , gradOnY(gradOnY_)
+        , gradOnZ(gradOnZ_)
+    {
+    }
+
+    double getSignedDis(const Eigen::Vector<double, dim>& point)
+    {
+        return disFunc(point);
+    }
+    double getGradOnX(const Eigen::Vector<double, dim>& point)
+    {
+        return gradOnX(point);
+    }
+    double getGradOnY(const Eigen::Vector<double, dim>& point)
+    {
+        return gradOnY(point);
+    }
+    double getGradOnZ(const Eigen::Vector<double, dim>& point)
+    {
+        return gradOnZ(point);
+    }
+
+private:
+    disCalculationPtr disFunc;
+    gradientDirXPtr gradOnX;
+    gradientDirXPtr gradOnY;
+    gradientDirXPtr gradOnZ;
+};
+
+}; // namespace ImplictSIM

include/Quadrature/Quadrature.hpp

@@ -0,0 +1,49 @@
+#pragma once
+#include "BgMesh/BgMesh.hpp"
+#include "Implict/ImplictModel.hpp"
+#include <algorithm>
+#include <Eigen/Eigen>
+#include <iostream>
+
+namespace ImplictSIM {
+
+enum class HeightDir { X = 0,
+    Y = 1,
+    Z = 2 };
+typedef std::pair<double, double> Interval;
+
+template <size_t dim>
+class VolumeIntegral {
+public:
+    typedef double (*integrandFunc)(const Eigen::Vector<double, dim>& point);
+    int dimension;
+    std::vector<ImplictModel<dim>> phiFuncs;
+    std::vector<int> si;
+    BgMesh<dim> meshU;
+    integrandFunc f;
+    bool flagS;
+    int GaussOrderQ;
+
+public:
+    VolumeIntegral(int dimension_, std::vector<ImplictModel<dim>>& phiFuncs_, const std::vector<int>& si_, const BgMesh<dim>& meshU_,
+        const integrandFunc f_, bool flagS_, int GaussOrderQ_)
+        : dimension(dimension_)
+        , phiFuncs(phiFuncs_)
+        , si(si_)
+        , meshU(meshU_)
+        , f(f_)
+        , flagS(flagS_)
+        , GaussOrderQ(GaussOrderQ_)
+    {
+    }
+
+    void findRoots(ImplictModel<dim> phiFunc,int si, HeightDir k, 
+    Interval interval, Eigen::Vector<double, dim > x, std::vector<double> roots,double eps);
+
+    void getOneIntegrationOnInterval(ImplictModel<dim> phiFunc, int si,
+        HeightDir k, Interval interval, Eigen::Vector<double, dim> x, std::vector<double>& hDirPoints);
+
+    void getIntegration();
+};
+
+}; // namespace ImplictSIM

src/Quadrature/Quadrature.cpp

@@ -0,0 +1,71 @@
+#include "Quadrature/Quadrature.hpp"
+
+namespace ImplictSIM {
+template <size_t dim>
+void VolumeIntegral<dim>::findRoots(ImplictModel<dim> phiFunc, int si, HeightDir k,
+    Interval interval, Eigen::Vector<double, dim> x, std::vector<double> roots, double eps)
+{
+    Eigen::Vector<double, dim> ek = Eigen::Vector<double, dim>::Zero();
+    if (k == HeightDir::X) {
+        ek(0) = 1.0;
+    } else if (K == HeightDir::Y) {
+        ek(1) = 1.0;
+    } else if (K == HeightDir::Z) {
+        ek(2) = 1.0;
+    }
+    // just binary search
+    double leftBound = interval.first, rightBound = interval.second;
+    double y;
+    roots.clear();
+    bool flag = false;
+    while (leftBound <= rightBound) {
+        y = (leftBound + rightBound) / 2.0;
+        if (phiFunc(x + y * ek) < -eps) {
+            leftBound = y;
+        } else if (phiFunc(x + y * ek) > eps) {
+            rightBound = y;
+        } else {
+            flag = true;
+            roots.push_back(y);
+            break;
+        }
+    }
+    roots.push_back(interval.first);
+    roots.push_back(interval.second);
+    sort(roots.begin(), root.end());
+    return;
+}
+
+template <size_t dim>
+void VolumeIntegral<dim>::getOneIntegrationOnInterval(ImplictModel<dim> phiFunc, int si,
+    HeightDir k, Interval interval, Eigen::Vector<double, dim> x, std::vector<double>& hDirPoints)
+{
+    Eigen::Vector<double, dim> ek = Eigen::Vector<double, dim>::Zero();
+    if (k == HeightDir::X) {
+        ek(0) = 1.0;
+    } else if (K == HeightDir::Y) {
+        ek(1) = 1.0;
+    } else if (K == HeightDir::Z) {
+        ek(2) = 1.0;
+    }
+    std::vector<double> roots;
+    findRoots(phiFunc, si, k, interval, x, root, 1e-8);
+    Eigen::Vector<double, dim> xc;
+    for (int i = 1)
+        for (int i = 1; i < roots.size(); i++) {
+            double len = roots[i] - root[i - 1];
+            xc = x + 0.5 * (root[i - 1] + root[i]) * ek;
+            if()
+        }
+}
+
+template <size_t dim>
+void VolumeIntegral<dim>::getIntegration()
+{
+    for (int curDim = dim; curDim >= 1; curDim--) {
+        if (curDim == 1) {
+        }
+    }
+}
+
+}; // namespace ImplictSIM

src/main.cpp

@@ -0,0 +1,7 @@
+#include <iostream>
+
+int main()
+{
+
+    return 0;
+}