Rename Context::Sequential to Context::Local

CMakeLists.txt

@@ -1,6 +1,7 @@
 cmake_minimum_required (VERSION 3.16)
 
 cmake_policy(SET CMP0054 NEW)
+cmake_policy(SET CMP0074 NEW)
 cmake_policy(SET CMP0076 NEW)
 
 # Prohibit in-source builds

examples/BoundaryOptimization/AcousticCloaking.cpp

@@ -29,14 +29,14 @@ const Attribute dBox = 16;
 const Attribute dSoundSoft = 11;
 const size_t maxIt = 2000;
 const Real epsilon = 1e-6;
-const Real ell = 1e-12;
+const Real ell = 1e-7;
 const Real tgv = 1e+12;
 const Real alpha = 4;
-const Real waveLength = 50;
+const Real waveLength = 20;
 const Real waveNumber = 2 * Math::Constants::pi() / waveLength;
 
 const VectorFunction xi = { 0, 0, 1 };
-const Real impedance = 1.0 / 10;
+const Real impedance = 1.0 / 20;
 const Real R = 50;
 
 using RealFES = P1<Real>;
@@ -48,13 +48,13 @@ using ShapeGradient = VectorGridFunction;
 template <class ProblemType>
 void solve(ProblemType& pb)
 {
-  Solver::UMFPack solver(pb);
-  solver.solve();
-  solver.printStatus();
-  // Solver::GMRES solver(pb);
-  // solver.setMaxIterations(2000);
-  // solver.setTolerance(1e-12);
+  // Solver::UMFPack solver(pb);
   // solver.solve();
+  // solver.printStatus();
+  Solver::GMRES solver(pb);
+  solver.setMaxIterations(3000);
+  solver.setTolerance(1e-12);
+  solver.solve();
 }
 
 int main(int, char**)
@@ -64,74 +64,67 @@ int main(int, char**)
   Threads::getGlobalThreadPool().reset(8);
   std::cout << Eigen::nbThreads() << std::endl;
 
+  MMG::Mesh miaow;
+  miaow.load("Scattered.mesh");
+
+
+  P1 fes(miaow);
+  GridFunction gf(fes);
+  gf.load("Scattered.gf");
+
+  GridFunction diff(fes);
+
+  diff = abs(gf - 1);
+  diff.save("Diff.gf");
+
+
+  std::exit(1);
+
   Alert::Info() << "Loading mesh..." << Alert::Raise;
   MMG::Mesh mesh;
   mesh.load("AcousticCloaking.medit.mesh", IO::FileFormat::MEDIT);
 
-  Real hmax = waveLength / 2;
-  Real hmin = waveLength / 32;
-  Real hausd = waveLength / 64;
-  Real hgrad = 1.2;
-
-  Alert::Info() << "Initializing sound-soft region..." << Alert::Raise;
-  {
-    P1 vh(mesh);
 
-    GridFunction dist(vh);
-    dist = [&](const Point& p)
-    {
-      Math::SpatialVector<Real> c(3);
-      c << -45.189708, 0, -11.330562;
-      return (p - c).norm() - 2;
-    };
+  Real hmax = waveLength / 3;
+  // Real hmin = 0.5;
+  // Real hausd = 0.1;
 
-    mesh.save("Dist.mesh", IO::FileFormat::MEDIT);
-    dist.save("Dist.sol", IO::FileFormat::MEDIT);
-
-    mesh = MMG::ImplicitDomainMesher().setAngleDetection(false)
-                                      .split(SoundSoft, { SoundSoft, SoundHard })
-                                      .split(SoundHard, { SoundSoft, SoundHard })
-                                      .noSplit(PML)
-                                      .noSplit(dBox)
-                                      .setHMax(hmax)
-                                      .setHMin(hmin)
-                                      .setGradation(hgrad)
-                                      .setHausdorff(hausd)
-                                      .surface()
-                                      .discretize(dist);
-  }
+  Real hmin = waveLength / 32;
+  Real hausd = waveLength / 64;
 
   mesh.save("Omega0.mfem.mesh", IO::FileFormat::MFEM);
   mesh.save("Omega0.mesh", IO::FileFormat::MEDIT);
 
   std::ofstream fObj("obj.txt");
+  std::ofstream fObjMax("obj_max.txt");
   size_t i = 0;
   size_t regionCount;
   while (i < maxIt)
   {
-    //bool topologicalStep = i < 5 || (i < 200  && i % 20 == 0);
-    bool topologicalStep = true;
+    bool topologicalStep = i < 6 || (i < 50  && i % 10 == 0);
+    // bool topologicalStep = false;
     bool geometricStep = !topologicalStep;
 
     const Real k = 0.5 * (hmax + hmin);
-    const Real dt = 2 * hmin;
-    const Real radius = 2 * hmin;
+    const Real dt = 4 * hmin;
+    const Real radius = 2;
     Alert::Info() << "Iteration: " << i                         << Alert::NewLine
                   << "HMax:      " << Alert::Notation(hmax)     << Alert::NewLine
                   << "HMin:      " << Alert::Notation(hmin)     << Alert::NewLine
                   << "Hausdorff: " << Alert::Notation(hausd)    << Alert::NewLine
-                  << "HGrad:     " << Alert::Notation(hgrad)    << Alert::NewLine
                   << "ell:     " << Alert::Notation(ell)    << Alert::NewLine
                   << "dt:        " << Alert::Notation(dt)       << Alert::Raise;
 
     Alert::Info() << "Optimizing the domain..." << Alert::Raise;
     MMG::Optimizer().setHMax(hmax)
                     .setHMin(hmin)
                     .setHausdorff(hausd)
-                    .setGradation(hgrad)
                     .setAngleDetection(false)
+                    .setGradation(1.2)
                     .optimize(mesh);
 
+    mesh.save("Optimized.mesh", IO::FileFormat::MEDIT);
+
     Alert::Info() << "Computing required connectivity..." << Alert::Raise;
     mesh.getConnectivity().compute(2, 3); // Computes boundary
     mesh.getConnectivity().compute(2, 2);
@@ -198,42 +191,52 @@ int main(int, char**)
     Alert::Info() << "Solving..." << Alert::Raise;
     solve(state);
 
+    Math::SpatialVector<Real> c{{-50, 0, 0}};
+
     GridFunction total(rfes);
-    total = [&](const Point& x) { return std::norm(wave(x) + u.getSolution()(x)); };
+    total = [&](const Point& x)
+    {
+      return std::norm(wave(x)+ u.getSolution()(x));
+    };
+
     total.save("Total.gf");
     rfes.getMesh().save("Total.mesh");
 
     GridFunction pressure(rfes);
-    Math::SpatialVector<Real> c({{-50, 0, 0}});
-    pressure = [&](const Point& p)
+    pressure = [&](const Point& x)
     {
-      if ((p - c).norm() < 75)
-        return std::norm(u.getSolution()(p));
-      else return 0.0;
+      return std::norm(u.getSolution()(x));
     };
     pressure.save("Scattered.gf");
     rfes.getMesh().save("Scattered.mesh");
 
     Alert::Info() << "Computing objective..." << Alert::Raise;
     pressure.setWeights();
-    const Real J = 0.5 * Integral(pressure).compute();
+    const Real J = 0.5 * Integral(pressure).compute() / mesh.getVolume();
     const Real area = mesh.getArea(SoundSoft);
     const Real objective = J + ell * area;
 
+    const Real objectiveMax = pressure.max() + ell * area;
+
     Alert::Info() << "Objective: " << Alert::Notation(objective) << Alert::NewLine
                   << "J: " << Alert::Notation(J) << Alert::NewLine
                   << "Area: " << Alert::Notation(area) << Alert::NewLine
+                  << "Objective Max: " << Alert::Notation(objectiveMax) << Alert::NewLine
                   << Alert::Raise;
+
     fObj << objective << "\n";
     fObj.flush();
 
+    fObjMax << objectiveMax << "\n";
+    fObjMax.flush();
+
     Alert::Info() << "Assembling adjoint equation..." << Alert::Raise;
     TrialFunction p(cfes);
     TestFunction  q(cfes);
     Problem adjoint(p, q);
     adjoint = Integral(Grad(p), Grad(q))
             - waveNumber * waveNumber * Integral(p, q).over(SoundSoft)
-            + Integral(u.getSolution(), q)
+            + Integral(u.getSolution() / mesh.getVolume(), q)
             + Complex(0, -waveNumber / impedance) * BoundaryIntegral(p, q).over(SoundSoft)
             + Complex(1 / R, -waveNumber) * FaceIntegral(p, q).over(PML)
             ;
@@ -256,14 +259,19 @@ int main(int, char**)
       Problem topo(s, t);
       topo = alpha * alpha * Integral(Grad(s), Grad(t))
            + Integral(s, t)
-           - Integral(
-               Im(waveNumber / impedance * Conjugate(u.getSolution() + wave) * p.getSolution()), t).over(SoundHard)
+           + Integral(
+               Im(waveNumber / impedance * Conjugate(u.getSolution() + wave) *
+                   p.getSolution()), t).over(SoundHard)
            + tgv * Integral(s, t).over(dBox, SoundSoft);
 
       topo.assemble();
 
       solve(topo);
 
+      GridFunction norm(drfes);
+      norm = Abs(s.getSolution());
+      s.getSolution() /= norm.max();
+
       s.getSolution().save("Topo.gf");
       drfes.getMesh().save("Topo.mesh");
 
@@ -296,6 +304,7 @@ int main(int, char**)
               - FaceIntegral(
                   Im(waveNumber / impedance * Conjugate(u.getSolution() + wave) * p.getSolution()),
                   Dot(conormal, w)).over(dSoundSoft)
+              + ell * FaceIntegral(conormal, w).over(dSoundSoft)
               ;
 
       hilbert.assemble();
@@ -327,9 +336,10 @@ int main(int, char**)
                                       .noSplit(dBox)
                                       .setHMax(hmax)
                                       .setHMin(hmin)
-                                      .setGradation(hgrad)
+                                      .setGradation(1.2)
                                       .setHausdorff(hausd)
                                       .surface()
+                                      .setBoundaryReference(10)
                                       .discretize(workaround);
 
     Alert::Info() << "Saving files..." << Alert::Raise;

examples/BoundaryOptimization/Tools.h

@@ -19,7 +19,7 @@ namespace Rodin::Examples::BoundaryOptimization
           Geometry::Polytope::getVertices(
             Geometry::Polytope::Type::Point).col(0), it->getCoordinates());
       const Real tp = tf(p);
-      if (tp > 0 && (tp / tc) > (1 - 1e-12))
+      if (tp > 1e-12 && (tp / tc) > (1 - 1e-12))
       {
         cs.emplace_back(std::move(p));
         break;

examples/BoundaryOptimization/WaterTank.cpp

@@ -116,6 +116,13 @@ int main(int, char**)
   Real hausd = 0.5 * hmin;
   Real hgrad = 1.2;
 
+  MMG::Optimizer().setHMax(hmax)
+                  .setHMin(hmin)
+                  .setHausdorff(hausd)
+                  .setGradation(hgrad)
+                  .setAngleDetection(true)
+                  .optimize(mesh);
+
   Alert::Info() << "Initializing unsupported region..." << Alert::Raise;
   {
     P1 vh(mesh);
@@ -128,7 +135,7 @@ int main(int, char**)
       return (p - c).norm() - 6.5;
     };
 
-    mesh = MMG::ImplicitDomainMesher().setAngleDetection(false)
+    mesh = MMG::ImplicitDomainMesher().setAngleDetection(true)
                                       .split(Gamma, { Unsupported, Gamma })
                                       .split(Unsupported, { Unsupported, Gamma })
                                       .setHMax(hmax)
@@ -375,8 +382,8 @@ int main(int, char**)
                                         .split(Gamma, { Support, Gamma })
                                         .noSplit(Unsupported)
                                         .setHMax(hmax)
-                                        // .setHMin(hmin)
-                                        // .setGradation(hgrad)
+                                        .setHMin(hmin)
+                                        .setGradation(hgrad)
                                         .setHausdorff(hausd)
                                         .surface()
                                         .discretize(workaround);

examples/Geometry/CMakeLists.txt

@@ -48,3 +48,8 @@ add_executable(UniformGrid UniformGrid.cpp)
 target_link_libraries(UniformGrid
   PUBLIC
   Rodin::Geometry)
+
+add_executable(UniformRefinement UniformRefinement.cpp)
+target_link_libraries(UniformRefinement
+  PUBLIC
+  Rodin::Geometry)

examples/Geometry/UniformGrid.cpp

@@ -14,7 +14,7 @@ int main(int, char**)
 {
   constexpr size_t n = 64;
   Mesh mesh;
-  mesh = SequentialMesh::UniformGrid(Polytope::Type::Triangle, { n, n });
+  mesh = LocalMesh::UniformGrid(Polytope::Type::Triangle, { n, n });
   mesh.save("UniformGrid.mesh");
   return 0;
 }

examples/Geometry/UniformRefinement.cpp

@@ -0,0 +1,20 @@
+/*
+ *          Copyright Carlos BRITO PACHECO 2021 - 2022.
+ * Distributed under the Boost Software License, Version 1.0.
+ *       (See accompanying file LICENSE or copy at
+ *          https://www.boost.org/LICENSE_1_0.txt)
+ */
+#include <Rodin/Geometry.h>
+#include <Rodin/Alert.h>
+
+using namespace Rodin;
+using namespace Geometry;
+
+int main(int, char**)
+{
+  Mesh mesh;
+  mesh = mesh.UniformGrid(Polytope::Type::Triangle, { 16, 16 });
+
+  return 0;
+}
+

examples/IO/MFEM.cpp

@@ -9,7 +9,7 @@ using namespace Rodin::Variational;
 int main(int, char** argv)
 {
   Mesh mesh =
-    Mesh<Rodin::Context::Sequential>::Builder()
+    Mesh<Rodin::Context::Local>::Builder()
     .initialize(2)
     .nodes(4)
     .vertex({0, 0})

examples/PDEs/Helmholtz.cpp

@@ -47,7 +47,7 @@ int main(int, char**)
   ComplexFunction f =
     [&](const Point& p)
     {
-      return Complex(p.x(), 1) * waveNumber * waveNumber * cos(waveNumber * p.x()
+      return Complex(1, 1) * waveNumber * waveNumber * cos(waveNumber * p.x()
           ) * cos(waveNumber * p.y());
     };
 

examples/PDEs/Poisson.cpp

@@ -4,6 +4,7 @@
  *       (See accompanying file LICENSE or copy at
  *          https://www.boost.org/LICENSE_1_0.txt)
  */
+#include <Rodin/Types.h>
 #include <Rodin/Solver.h>
 #include <Rodin/Geometry.h>
 #include <Rodin/Variational.h>
@@ -17,21 +18,27 @@ int main(int, char**)
   // Build a mesh
   Mesh mesh;
   mesh = mesh.UniformGrid(Polytope::Type::Triangle, { 16, 16 });
-  mesh.getConnectivity().compute(1, 2);
 
   // Functions
-  P1 vh(mesh);
+  P1<Complex> vh(mesh);
 
-  TrialFunction u(vh);
-  TestFunction  v(vh);
+  GridFunction miaow(vh);
+  miaow = [](const Point& p){ return Complex(0, 1) * Complex(p.x(), 2 * p.y()); };
+  miaow.setWeights();
 
-  // Define problem
-  Problem poisson(u, v);
-  poisson = Integral(Grad(u), Grad(v))
-          - Integral(v)
-          + DirichletBC(u, Zero());
+  Complex s = Integral(miaow).compute();
+  std::cout << s << std::endl;
 
-  poisson.assemble();
+  // TrialFunction u(vh);
+  // TestFunction  v(vh);
+
+  // // Define problem
+  // Problem poisson(u, v);
+  // poisson = Integral(Grad(u), Grad(v))
+  //         - Integral(v)
+  //         + DirichletBC(u, Zero());
+
+  // poisson.assemble();
   // Solver::CG solver;
   // poisson.solve(solver);