Example for Geometric Advection, Boolean Operations, and Iterative Advection in ViennaLS

iterativeAfterGeometricAdvect.cpp

#include <array>

#include <lsAdvect.hpp>
#include <lsBooleanOperation.hpp>
#include <lsDomain.hpp>
#include <lsGeometricAdvect.hpp>
#include <lsGeometricAdvectDistributions.hpp>
#include <lsGeometries.hpp>
#include <lsMakeGeometry.hpp>
#include <lsMesh.hpp>
#include <lsPrune.hpp>
#include <lsToMesh.hpp>
#include <lsToSurfaceMesh.hpp>
#include <lsVTKWriter.hpp>

template <class T, int D>
static void outputLS(lsSmartPointer<lsDomain<T, D>> ls, std::string fileName) {
  auto mesh = lsSmartPointer<lsMesh<T>>::New();
  lsToSurfaceMesh<T, D>(ls, mesh).apply();
  lsVTKWriter<T>(mesh, fileName).apply();
  // outputting the LS values as a mesh is often useful for debugging
  // lsToMesh<T, D>(ls, mesh, false).apply();
  // std::string newFileName = fileName.substr(0, fileName.size() - 4) +
  //                           "_points" + fileName.substr(fileName.size() - 4);
  // lsVTKWriter<T>(mesh, newFileName).apply();
}

template <int Dimensions>
static void
BooleanAddition(lsSmartPointer<lsDomain<double, Dimensions>> &target,
                lsSmartPointer<lsDomain<double, Dimensions>> &source) {
  lsBooleanOperation<double, Dimensions>(target, source,
                                           lsBooleanOperationEnum::UNION).apply();
  // do not need to prune as this is done inside BooleanOps already
  // lsPrune<double, Dimensions> thin_out(target);
  // thin_out.apply();
}

template <int Dimensions>
static void
BooleanSubtraction(lsSmartPointer<lsDomain<double, Dimensions>> &target,
                   lsSmartPointer<lsDomain<double, Dimensions>> &source) {
  lsBooleanOperation<double, Dimensions>(
      target, source, lsBooleanOperationEnum::RELATIVE_COMPLEMENT).apply();
}

template <class T> class VelocityField : public lsVelocityField<T> {
  T getScalarVelocity(const std::array<T, 3> &coordinate, int material,
                      const std::array<T, 3> &normalVector,
                      unsigned long pointId) final {
    return 0.1;
  }
};

int main() {
  omp_set_num_threads(1);

  constexpr int D = 2;
  typedef double NumericType;
  double gridDelta = 0.5;

  // Set up domain properties
  lsSmartPointer<lsDomain<NumericType, D>> substrate;
  {
    double extent = 4.5;
    double bounds[2 * D] = {-extent, extent, -extent, extent};

    typename lsDomain<NumericType, D>::BoundaryType boundaryCons[D];
    boundaryCons[0] =
        lsDomain<NumericType, D>::BoundaryType::REFLECTIVE_BOUNDARY;
    boundaryCons[1] = lsDomain<NumericType, D>::BoundaryType::INFINITE_BOUNDARY;

    if constexpr (D == 3) {
      bounds[4] = -extent;
      bounds[5] = extent;
      boundaryCons[1] =
          lsDomain<NumericType, D>::BoundaryType::REFLECTIVE_BOUNDARY;
      boundaryCons[2] =
          lsDomain<NumericType, D>::BoundaryType::INFINITE_BOUNDARY;
    }

    // Set up level set with these values, grid should just be used to create new
    // LS on the same grid
    substrate = lsSmartPointer<lsDomain<NumericType, D>>::New(
        bounds, boundaryCons, gridDelta);
  }

  // now create mask LS
  // just use copy of other grid to initialise empty mask LS
  auto maskLS = lsSmartPointer<lsDomain<NumericType, D>>::New(substrate->getGrid());
  {
    // mask should extend slightly into substrate at this low resolution
    // due to numerical issues
    std::array<NumericType, D> min = {-1., 19.5};
    std::array<NumericType, D> max = {1., 27.5};
    lsMakeGeometry(maskLS, lsSmartPointer<lsBox<NumericType, D>>::New(min.data(), max.data())).apply();
  }

  // now set up substrate
  {
    std::array<double, D> p = {0.0, 20.0};
    std::array<double, D> n = {0.0, 1.0};
    lsMakeGeometry<double, D>(
        substrate, lsSmartPointer<lsPlane<double, D>>::New(p.data(), n.data()))
        .apply();
    BooleanAddition(substrate, maskLS);

    // now print LSs as they are initially
    outputLS(maskLS, "mask_initial.vtp");
    outputLS(substrate, "subs_initial.vtp");
  }

  // do the geometric advection
  {
    auto dist = lsSmartPointer<lsBoxDistribution<double, D>>::New(
      std::array<double, 3>({-gridDelta, -5.0, -gridDelta}), gridDelta);
    lsGeometricAdvect<NumericType, D> ga;
    ga.setMaskLevelSet(maskLS);
    ga.setLevelSet(substrate);
    ga.setAdvectionDistribution(dist);
    ga.apply();

    outputLS(maskLS, "mask_etched.vtp");
    outputLS(substrate, "subs_etched.vtp");
  }

  // subtact mask from substrate
  lsBooleanOperation(substrate, maskLS, lsBooleanOperationEnum::RELATIVE_COMPLEMENT).apply();
  // BooleanSubtraction(substrate, maskLS);
  outputLS(substrate, "subs_booled.vtp");

  // now deposit isotropically with iterative advection
  {
    auto velocityField = lsSmartPointer<VelocityField<double>>::New();

    lsAdvect<double, D> advect;
    // EO_1st is already default, so no need to set it explicitly
    // advect.setIntegrationScheme(
    //     lsIntegrationSchemeEnum::ENGQUIST_OSHER_1ST_ORDER);
    advect.insertNextLevelSet(substrate);
    advect.setTimeStepRatio(0.1);
    advect.setAdvectionTime(5.0);
    advect.setVelocityField(velocityField);
    advect.apply();
    outputLS(substrate, "subs_deposited.vtp");
  }

  return 0;
}