Develop

examples/BoundaryOptimization/BoundaryHeatConductivity.cpp

@@ -31,8 +31,8 @@ static constexpr Scalar ell = 1;
 static constexpr Scalar radius = 0.02;
 static constexpr Scalar tgv = std::numeric_limits<float>::max();
 
-using ScalarFES = P1<Scalar, Context::Serial>;
-using VectorFES = P1<Math::Vector, Context::Serial>;
+using ScalarFES = P1<Scalar, Context::Sequential>;
+using VectorFES = P1<Math::Vector, Context::Sequential>;
 using ScalarGridFunction = GridFunction<ScalarFES>;
 using VectorGridFunction = GridFunction<VectorFES>;
 using ShapeGradient = VectorGridFunction;

examples/Geometry/UniformGrid.cpp

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

examples/IO/MFEM.cpp

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

examples/IntegralEquations/EquilibriumDistribution.cpp

@@ -15,44 +15,90 @@ using namespace Rodin::Variational;
 inline
 Scalar K(const Point& x, const Point& y)
 {
-  return 1. / (4 * M_PI * (x - y).norm());
+  return 1. / (4 * M_PI * ((x - y).norm()));
+}
+
+inline
+Scalar one(const Point& x, const Point& y)
+{
+  return 1;
+}
+// inline
+// Scalar K(const Point& x, const Point& y)
+// {
+//   const Scalar n = (x - y).norm();
+//   return 1. / (4 * M_PI * (n * n * n));
+// }
+
+inline
+Scalar Ke(const Point& x, const Point& y)
+{
+  return 1. / (4 * M_PI * ((x - y).norm()));
 }
 
 inline
 Scalar exact(const Point& x)
 {
-  return 4. / (M_PI * std::sqrt(1 - x.squaredNorm()));
+  if (abs(1 - x.squaredNorm()) < 0.001)
+    return 4. / (M_PI * std::sqrt(abs(0.001)));
+  else
+    return 4. / (M_PI * std::sqrt(abs(1 - x.squaredNorm())));
 }
 
 int main(int, char**)
 {
   Mesh mesh;
   mesh.load("D1.mesh");
-  // mesh.save("miaow.medit.mesh", IO::FileFormat::MEDIT);
+  // mesh.load("miaow.medit.o.mesh", IO::FileFormat::MEDIT);
   mesh.getConnectivity().compute(1, 2);
 
   P1 fes(mesh);
-
   TrialFunction u(fes);
   TestFunction  v(fes);
+
   DenseProblem eq(u, v);
-  eq = Integral(Potential(K, u), v)
-     - Integral(v);
+  eq = Integral(0.00001 * Grad(u), Grad(v))
+     + Integral(Potential(K, u), v)
+     - Integral(v)
+     ;
 
   std::cout << "assemblage\n";
   eq.assemble();
 
+  std::cout << "save\n";
+  std::ofstream file("matrix.csv");
+  file << eq.getStiffnessOperator().format(
+      Eigen::IOFormat(Eigen::StreamPrecision, 0, ", ", "\n"));
+  file.close();
+
   std::cout << "resolution\n";
   Solver::LDLT solver;
   eq.solve(solver);
 
   u.getSolution().save("u.gf");
   mesh.save("u.mesh");
 
-  // GridFunction phi(fes);
-  // phi = [](const Point& p) { return 4. / (M_PI * std::sqrt(1 - p.squaredNorm())); };
-  // phi.projectOnBoundary(ScalarFunction(0));
-  // phi.save("phi.gf");
+  std::cout << "average:\n";
+  u.getSolution().setWeights();
+  std::cout << Integral(u.getSolution()) << std::endl;
+
+  GridFunction one(fes);
+  one = ScalarFunction(1);
+
+  std::cout << "exact\n";
+  GridFunction ex(fes);
+  ex = exact;
+  ex.save("exact.gf");
+
+  std::cout << "phi\n";
+  GridFunction phi(fes);
+  phi = Potential(K, u.getSolution());
+  phi.save("phi.gf");
+
+  std::cout << "potential\n";
+  phi.setWeights();
+  ex.setWeights();
+  std::cout << Integral(phi).compute() << std::endl;
 
   return 0;
 }

examples/QF/GrundmannMoller.cpp

@@ -12,15 +12,25 @@ using namespace Rodin::Geometry;
 
 int main(int, char**)
 {
-  QF::GrundmannMoller qf(1, Polytope::Type::Triangle);
-  std::cout << "Order: " << qf.getOrder() << std::endl;
-  std::cout << "Size: " << qf.getSize() << std::endl;
-  for (size_t i = 0; i < qf.getSize(); i++)
+  QF::GrundmannMoller qfe(1, Polytope::Type::Segment);
+  std::cout << "Order: " << qfe.getOrder() << std::endl;
+  std::cout << "Size: " << qfe.getSize() << std::endl;
+  for (size_t i = 0; i < qfe.getSize(); i++)
   {
     std::cout << "Weight:\n";
-    std::cout << qf.getWeight(i) << std::endl;
-    std::cout << "Point:\n" << qf.getPoint(i) << std::endl;
+    std::cout << qfe.getWeight(i) << std::endl;
+    std::cout << "Point:\n" << qfe.getPoint(i) << std::endl;
   }
 
+  // QF::GrundmannMoller qf(1, Polytope::Type::Triangle);
+  // std::cout << "Order: " << qf.getOrder() << std::endl;
+  // std::cout << "Size: " << qf.getSize() << std::endl;
+  // for (size_t i = 0; i < qf.getSize(); i++)
+  // {
+  //   std::cout << "Weight:\n";
+  //   std::cout << qf.getWeight(i) << std::endl;
+  //   std::cout << "Point:\n" << qf.getPoint(i) << std::endl;
+  // }
+
   return 0;
 }

examples/ShapeOptimization/CMakeLists.txt

@@ -16,15 +16,15 @@ target_link_libraries(LevelSetCantilever2D
   Rodin::External::MMG
   )
 
-# add_executable(LevelSetCantilever3D LevelSetCantilever3D.cpp)
-# target_link_libraries(LevelSetCantilever3D
-#   PUBLIC
-#   Rodin::Geometry
-#   Rodin::Solver
-#   Rodin::Variational
-#   Rodin::External::MMG
-#   )
-# 
+add_executable(LevelSetCantilever3D LevelSetCantilever3D.cpp)
+target_link_libraries(LevelSetCantilever3D
+  PUBLIC
+  Rodin::Geometry
+  Rodin::Solver
+  Rodin::Variational
+  Rodin::External::MMG
+  )
+
 # add_executable(LevelSetMast2D LevelSetMast2D.cpp)
 # target_link_libraries(LevelSetMast2D
 #   PUBLIC

examples/ShapeOptimization/LevelSetCantilever2D.cpp

@@ -15,7 +15,7 @@ using namespace Rodin::Geometry;
 using namespace Rodin::External;
 using namespace Rodin::Variational;
 
-using FES = VectorP1<Context::Serial>;
+using FES = VectorP1<Context::Sequential>;
 
 // Define interior and exterior for level set discretization
 static constexpr Attribute Interior = 1, Exterior = 2;
@@ -108,7 +108,7 @@ int main(int, char**)
     hilbert = Integral(alpha * Jacobian(g), Jacobian(w))
             + Integral(g, w)
             - FaceIntegral(Dot(Ae, e) - ell, Dot(n, w)).over(Gamma)
-            + DirichletBC(g, VectorFunction{0, 0}).on(GammaN);
+            + DirichletBC(g, VectorFunction{0, 0, 0}).on(GammaN);
     hilbert.solve(solver);
     const auto& dJ = g.getSolution();
 
@@ -136,7 +136,7 @@ int main(int, char**)
 
     th = MMG::ImplicitDomainMesher().split(Interior, {Interior, Exterior})
                                     .split(Exterior, {Interior, Exterior})
-                                    .setRMC(1e-3)
+                                    .setRMC(1e-6)
                                     .setHMax(hmax)
                                     .setAngleDetection(false)
                                     .setBoundaryReference(Gamma)