Doc

doc/Namespaces.dox

@@ -191,6 +191,47 @@
   where @f$ \mathbf{A} @f$ represents the stiffness operator.
  */
 
+/**
+  @dir Rodin/Models
+  @brief Namespace @ref Rodin::Models
+ */
+
+/**
+  @namespace Rodin::Models
+  @brief Module which provides implementation of various variational models.
+ */
+
+/**
+  @dir Rodin/Models/Distance
+  @brief Namespace @ref Rodin::Models::Distance
+ */
+
+/**
+  @namespace Rodin::Models::Distance
+  @brief Module which provides models for computation of the distance function.
+ */
+
+/**
+  @dir Rodin/Models/ShapeOptimization
+  @brief Namespace @ref Rodin::Models::ShapeOptimization
+ */
+
+/**
+  @namespace Rodin::Models::ShapeOptimization
+  @brief Module which provides models and tools for geometrical shape optimization.
+ */
+
+/**
+  @dir Rodin/Models/Advection
+  @brief Namespace @ref Rodin::Models::Advection
+ */
+
+/**
+  @namespace Rodin::Models::Advection
+  @brief Module which provides models for the resolution of unsteady advection
+  equations.
+ */
+
 /**
   @dir Rodin/QF
   @brief Namespace @ref Rodin::QF

doc/References.bib

@@ -73,3 +73,33 @@ @article{logg2009efficient
   year={2009},
   publisher={Inderscience Publishers}
 }
+
+@inproceedings{belyaev2015variational,
+  title={On variational and PDE-based distance function approximations},
+  author={Belyaev, Alexander G and Fayolle, Pierre-Alain},
+  booktitle={Computer Graphics Forum},
+  volume={34},
+  number={8},
+  pages={104--118},
+  year={2015},
+  organization={Wiley Online Library}
+}
+
+@article{tucker1998assessment,
+  title={Assessment of geometric multilevel convergence robustness and a wall distance method for flows with multiple internal boundaries},
+  author={Tucker, PG},
+  journal={Applied Mathematical Modelling},
+  volume={22},
+  number={4-5},
+  pages={293--311},
+  year={1998},
+  publisher={Elsevier}
+}
+
+@inproceedings{spalding1994calculation,
+  title={Calculation of turbulent heat transfer in cluttered spaces},
+  author={Spalding, DB},
+  booktitle={Proceedings of the 10th International Heat Transfer Conference, Brighton, UK},
+  pages={14--18},
+  year={1994}
+}

examples/ShapeOptimization/LevelSetCantilever2D.cpp

@@ -73,7 +73,6 @@ int main(int, char**)
 
     Alert::Info() << "   | Trimming mesh." << Alert::Raise;
     SubMesh trimmed = th.trim(Exterior);
-    trimmed.getConnectivity().compute(1, 2);
 
     Alert::Info() << "   | Building finite element spaces." << Alert::Raise;
     const size_t d = 2;

examples/SurfaceEvolution/ConormalAdvection/;

@@ -0,0 +1,270 @@
+#include <thread>
+#include <queue>
+#include <Rodin/Math.h>
+#include <Rodin/Solver.h>
+#include <Rodin/Geometry.h>
+#include <Rodin/Variational.h>
+#include <RodinExternal/MMG.h>
+
+using namespace std;
+using namespace Rodin;
+using namespace Rodin::Geometry;
+using namespace Rodin::Variational;
+using namespace Rodin::External;
+
+constexpr Geometry::Attribute sphereCap = 3;
+constexpr char meshFile[] =
+  "../resources/examples/SurfaceEvolution/FirePropagation/GroundPlane.mfem.mesh";
+
+class ThreadPool {
+public:
+    ThreadPool(size_t numThreads) : stop(false) {
+        for (size_t i = 0; i < numThreads; ++i) {
+            workers.emplace_back([this] {
+                while (true) {
+                    std::function<void()> task;
+
+                    {
+                        std::unique_lock<std::mutex> lock(queueMutex);
+                        condition.wait(lock, [this] { return stop || !tasks.empty(); });
+
+                        if (stop && tasks.empty()) {
+                            return;
+                        }
+
+                        task = std::move(tasks.front());
+                        tasks.pop();
+                    }
+
+                    task();
+                }
+            });
+        }
+    }
+
+    template <typename Func>
+    void enqueue(Func func) {
+        {
+            std::unique_lock<std::mutex> lock(queueMutex);
+            tasks.emplace(func);
+        }
+        condition.notify_one();
+    }
+
+    ~ThreadPool() {
+        {
+            std::unique_lock<std::mutex> lock(queueMutex);
+            stop = true;
+        }
+        condition.notify_all();
+        for (std::thread &worker : workers) {
+            worker.join();
+        }
+    }
+
+private:
+    std::vector<std::thread> workers;
+    std::queue<std::function<void()>> tasks;
+
+    std::mutex queueMutex;
+    std::condition_variable condition;
+    bool stop;
+};
+
+struct Experiment
+{
+  const double azimuth;
+  const double hmax;
+  const double c;
+  const double T;
+  const double hausd;
+};
+
+std::vector<Experiment> experiments = {};
+
+void run(size_t expId, const std::vector<Experiment>& experiments);
+
+double phi(double t, const Point& p, const double azimuth);
+
+int main(int argc, char** argv)
+{
+  size_t N = 32;
+  for (size_t i = 0; i < N; i++)
+    experiments.push_back(
+        {0.1, 0.02 + (1 - 0.02) * float(i) / float(N), 0.1, M_PI / 2 - 0.1, 0.01});
+
+  for (size_t i = 0; i < N; i++)
+    experiments.push_back(
+        {0.1, 0.02 + (1 - 0.02) * float(i) / float(N), 0.2, M_PI / 2 - 0.1, 0.01});
+
+  for (size_t i = 0; i < N; i++)
+    experiments.push_back(
+        {0.1, 0.02 + (1 - 0.02) * float(i) / float(N), 0.3, M_PI / 2 - 0.1, 0.01});
+
+  for (size_t i = 0; i < N; i++)
+    experiments.push_back(
+        {0.1, 0.02 + (1 - 0.02) * float(i) / float(N), 0.4, M_PI / 2 - 0.1, 0.01});
+
+  for (size_t i = 0; i < N; i++)
+    experiments.push_back(
+        {0.1, 0.02 + (1 - 0.02) * float(i) / float(N), 0.5, M_PI / 2 - 0.1, 0.01});
+
+  for (size_t i = 0; i < N; i++)
+    experiments.push_back(
+        {0.1, 0.02 + (1 - 0.02) * float(i) / float(N), 0.6, M_PI / 2 - 0.1, 0.01});
+
+  for (size_t i = 0; i < N; i++)
+    experiments.push_back(
+        {0.1, 0.02 + (1 - 0.02) * float(i) / float(N), 0.7, M_PI / 2 - 0.1, 0.01});
+
+  for (size_t i = 0; i < N; i++)
+    experiments.push_back(
+        {0.1, 0.02 + (1 - 0.02) * float(i) / float(N), 0.8, M_PI / 2 - 0.1, 0.01});
+
+  for (size_t i = 0; i < N; i++)
+    experiments.push_back(
+        {0.1, 0.02 + (1 - 0.02) * float(i) / float(N), 0.9, M_PI / 2 - 0.1, 0.01});
+
+  for (size_t i = 0; i < N; i++)
+    experiments.push_back(
+        {0.1, 0.02 + (1 - 0.02) * float(i) / float(N), 1.0, M_PI / 2 - 0.1, 0.01});
+
+  const size_t hwc = std::thread::hardware_concurrency();
+  const size_t n = hwc - 2;
+  ThreadPool threadPool(n);
+  for (size_t i = 0; i < experiments.size(); i++)
+    threadPool.enqueue([i]{ run(i, experiments); });
+
+  return 0;
+}
+
+void run(size_t experimentId, const std::vector<Experiment>& experiments)
+{
+  Experiment experiment = experiments[experimentId];
+
+  const bool physical = true;
+  double dt = NAN;
+  if (physical)
+    dt = experiment.c * experiment.hmax;
+  else
+    dt = experiment.c;
+
+  // Load mesh
+  MMG::Mesh th;
+  th.load(meshFile);
+
+  {
+    P1 vh(th);
+    GridFunction dist(vh);
+    dist =
+      [](const Point& p)
+      {
+        Math::SpatialVector c{{2.5e4, 2.5e4, 0}};
+        const Scalar d = (p - c).norm() - 0.5e4;
+        return d;
+      };
+    dist.save("dist.gf");
+    th.save("dist.mesh");
+    th = MMG::ImplicitDomainMesher().setAngleDetection(false)
+                                    .setHausdorff(100)
+                                    .discretize(dist);
+    th.save("implicit.mesh", IO::FileFormat::MEDIT);
+    std::exit(1);
+  }
+
+  std::string filename;
+  if (physical)
+    filename = "physical/L2ErrorPhysical_" + std::to_string(experimentId) + ".csv";
+  else
+    filename = "L2Error_" + std::to_string(experimentId) + ".csv";
+
+  std::ofstream fout(filename);
+  fout << "t,$\\mathcal{E}(t)$\n" << std::flush;
+  double t = 0;
+  size_t i = 0;
+  while (true)
+  {
+    // Build finite element space on the mesh
+    P1 vh(th);
+    P1 uh(th, th.getSpaceDimension());
+
+    // Distance the subdomain
+    GridFunction dist(vh);
+
+    if (i == 0)
+    {
+      dist = [&](const Point& p) { return phi(0, p, experiment.azimuth); };
+    }
+    else
+    {
+      dist = MMG::Distancer(vh).setInteriorDomain(sphereCap)
+                               .distance(th);
+
+      // Compute gradient of signed distance function
+      GridFunction conormal(uh);
+      auto gd = Grad(dist);
+      conormal = gd / Frobenius(gd);
+
+      // Advect
+      try
+      {
+        MMG::Advect(dist, conormal).step(std::min(dt, experiment.T - t));
+      }
+      catch (Alert::Exception& e)
+      {
+        Alert::Warning()
+          << "Advection failed for experiment " << experimentId
+          << Alert::NewLine
+          << "Retrying..."
+          << Alert::Raise;
+        continue;
+      }
+
+      t += std::min(dt, experiment.T - t);
+    }
+
+    auto phit = ScalarFunction(
+        [&](const Point& p) { return phi(t, p, experiment.azimuth); });
+
+    // Compute L2 error
+    GridFunction diff(vh);
+    diff = Pow(dist - phit, 2);
+    diff.setWeights();
+    double error = Integral(diff).compute();
+    fout << t << "," << error << '\n' << std::flush;
+
+    // Generate mesh to subdomain
+    th = MMG::ImplicitDomainMesher().setAngleDetection(false)
+                                    .setHMax(experiment.hmax)
+                                    .setHausdorff(experiment.hausd)
+                                    .discretize(dist);
+
+    MMG::Optimizer().setAngleDetection(false)
+                   .setHMax(experiment.hmax)
+                   .setHausdorff(experiment.hausd)
+                   .optimize(th);
+
+    // th.save("out/SphereCap.mfem." + std::to_string(i) + ".mesh");
+
+    if (t + std::numeric_limits<double>::epsilon() > experiment.T)
+      break;
+
+    // Alert::Info() << "l2: " << error << '\n' << Alert::Raise;
+    i++;
+  }
+}
+
+double phi(double t, const Point& p, const double azimuth)
+{
+  const double x = p.x();
+  const double y = p.y();
+  const double z = p.z();
+  assert(x != 0 && y != 0);
+  assert(x != 0 && y != 0 && z != 0);
+  const double alpha = std::acos(z);
+  assert(0 <= alpha && alpha <= M_PI);
+  const double beta = std::atan2(y, x) + M_PI;
+  assert(0 <= beta && beta <= 2 * M_PI);
+  return alpha - t - azimuth;
+}
+

examples/SurfaceEvolution/ConormalAdvection/CMakeLists.txt

@@ -15,3 +15,12 @@ target_link_libraries(Test
   Rodin::Variational
   Rodin::External::MMG
   )
+
+add_executable(Plane Plane.cpp)
+target_link_libraries(Plane
+  PUBLIC
+  Rodin::Solver
+  Rodin::Geometry
+  Rodin::Variational
+  Rodin::External::MMG
+  )

examples/SurfaceEvolution/ConormalAdvection/Main.cpp

@@ -17,7 +17,7 @@ constexpr Scalar hmin = 0.1 * hmax;
 constexpr size_t maxIt = 100;
 constexpr Geometry::Attribute sphereCap = 3;
 constexpr char meshFile[] =
-  "/Users/carlos/Projects/rodin/resources/examples/SurfaceEvolution/ConormalAdvection/SphereCap.medit.mesh";
+  "../resources/examples/SurfaceEvolution/ConormalAdvection/SphereCap.medit.mesh";
 
 int main()
 {