FIN DE ANNEE

CMakeLists.txt

@@ -91,6 +91,10 @@ if (RODIN_BUILD_SRC)
   find_package(Threads REQUIRED)
   add_subdirectory(third-party/BS_thread_pool)
 
+  if (RODIN_USE_OPENMP)
+    find_package(OpenMP REQUIRED)
+  endif()
+
   # ---- Boost ----
   # set(Boost_USE_STATIC_LIBS        OFF)
   # set(Boost_USE_STATIC_RUNTIME     OFF)

examples/IntegralEquations/CMakeLists.txt

@@ -1,5 +1,13 @@
-add_executable(EquilibriumDistribution EquilibriumDistribution.cpp)
-target_link_libraries(EquilibriumDistribution
+add_executable(NewtonianPotential NewtonianPotential.cpp)
+target_link_libraries(NewtonianPotential
+  PUBLIC
+  Rodin::Geometry
+  Rodin::Solver
+  Rodin::Variational
+  )
+
+add_executable(ElasticDistribution ElasticDistribution.cpp)
+target_link_libraries(ElasticDistribution
   PUBLIC
   Rodin::Geometry
   Rodin::Solver

examples/IntegralEquations/ElasticDistribution.cpp

@@ -0,0 +1,88 @@
+/*
+ *          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/Solver.h>
+#include <Rodin/Geometry.h>
+#include <Rodin/Variational.h>
+
+using namespace Rodin;
+using namespace Rodin::Geometry;
+using namespace Rodin::Variational;
+
+const Scalar lambda = 0.5769, mu = 0.3846;
+const Scalar nu = lambda / (2 * (lambda + mu));
+
+inline
+void K(Math::Matrix& res, const Point& x, const Point& y)
+{
+  const auto norm = (x - y).norm();
+  res.resize(3, 3);
+  const Scalar L00 = (1 - nu) / (2 * M_PI * mu * norm) + nu * (x(0) - y(0)) * (x(0) - y(0)) / (2 * M_PI * mu * norm * norm * norm);
+  const Scalar L01 = nu * (x(0) - y(0)) * (x(1) - y(1)) / (2 * M_PI * mu * norm * norm * norm);
+  const Scalar L02 = nu * (x(0) - y(0)) * (x(2) - y(2)) / (2 * M_PI * mu * norm * norm * norm);
+
+  const Scalar L10 = nu * (x(1) - y(1)) * (x(0) - y(0)) / (2 * M_PI * mu * norm * norm * norm);
+  const Scalar L11 = (1 - nu) / (2 * M_PI * mu * norm) + nu * (x(1) - y(1)) * (x(1) - y(1)) / (2 * M_PI * mu * norm * norm * norm);
+  const Scalar L12 = nu * (x(1) - y(1)) * (x(2) - y(2)) / (2 * M_PI * mu * norm * norm * norm);
+
+  const Scalar L20 = (1 - 2 * nu) * (x(0) - y(0)) / (4 * M_PI * mu * norm * norm);
+  const Scalar L21 = (1 - 2 * nu) * (x(1) - y(1)) / (4 * M_PI * mu * norm * norm);
+  const Scalar L22 = (1 - nu) / (2 * M_PI * mu * norm);
+
+  res << L00, L01, L02,
+         L10, L11, L12,
+         L20, L21, L22;
+}
+
+int main(int, char**)
+{
+  Mesh mesh;
+  mesh.load("../build/D1.mesh");
+  mesh.getConnectivity().compute(1, 2);
+
+  P1 fes(mesh, 3);
+
+  VectorFunction e1{0, 0, 1};
+
+  TrialFunction u(fes);
+  TestFunction  v(fes);
+  DenseProblem eq(u, v);
+  eq = 0.001 * LinearElasticityIntegral(u, v)(lambda, mu)
+     + Integral(Potential(K, u), v)
+     - Integral(e1, 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::CG<Math::Matrix, Math::Vector> cg;
+  eq.solve(cg);
+
+  u.getSolution().save("u.gf");
+  mesh.save("u.mesh");
+
+  // std::cout << "average:\n";
+  // u.getSolution().setWeights();
+  // std::cout << Integral(u.getSolution()) << std::endl;
+
+  std::cout << "phi\n";
+  GridFunction phi(fes);
+  phi = Potential(K, u.getSolution());
+  phi.save("phi.gf");
+
+  // std::cout << "potential\n";
+  // phi.setWeights();
+  // std::cout << Integral(phi).compute() / M_PI << std::endl;
+
+  return 0;
+}
+

examples/IntegralEquations/EquilibriumDistribution.cpp → examples/IntegralEquations/NewtonianPotential.cpp

@@ -67,20 +67,15 @@ int main(int, char**)
   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;
+  one.setWeights();
+  std::cout << Integral(phi).compute() / Integral(one) << std::endl;
 
   return 0;
 }

examples/PDEs/Elasticity.cpp

@@ -30,7 +30,7 @@ int main(int argc, char** argv)
   P1 fes(mesh, d);
 
   // Lamé coefficients
-  Scalar lambda = 0.5769, mu = 0.3846;
+  const Scalar lambda = 0.5769, mu = 0.3846;
 
   // Pull force
   VectorFunction f{0, -1};

examples/ShapeOptimization/LevelSetCantilever2D.cpp

@@ -61,10 +61,10 @@ int main(int, char**)
   Alert::Info() << "Saved initial mesh to Omega0.mesh" << Alert::Raise;
 
   // Solver
-  Solver::SparseLU solver;
+  Solver::CG solver;
 
   // Optimization loop
-  std::vector<double> obj;
+ std::vector<double> obj;
   std::ofstream fObj("obj.txt");
   for (size_t i = 0; i < maxIt; i++)
   {

examples/ShapeOptimization/LevelSetCantilever3D.cpp

@@ -35,7 +35,7 @@ static double hmax = 0.1;
 static double hmin = 0.1 * hmax;
 static double hausd = 0.5 * hmin;
 static size_t hmaxIt = maxIt / 2;
-const Scalar k = 4;
+const Scalar k = 1;
 const Scalar dt = k * (hmax - hmin);
 static double alpha = dt;
 
@@ -70,7 +70,7 @@ int main(int, char**)
   Alert::Info() << "Saved initial mesh to Omega0.mesh" << Alert::Raise;
 
   // Solver
-  Solver::SparseLU solver;
+  Solver::CG solver;
 
   // Optimization loop
   std::vector<double> obj;
@@ -120,7 +120,7 @@ int main(int, char**)
     TrialFunction g(vh);
     TestFunction  w(vh);
     Problem hilbert(g, w);
-    hilbert = Integral(alpha * alpha * Jacobian(g), Jacobian(w))
+    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, 0}).on(GammaN);
@@ -151,15 +151,13 @@ int main(int, char**)
                                     .setRMC(1e-5)
                                     .setHMax(hmax)
                                     .setHMin(hmin)
-                                    .setHausdorff(hausd)
                                     .setAngleDetection(false)
                                     .setBoundaryReference(Gamma)
                                     .setBaseReferences(GammaD)
                                     .discretize(dist);
 
     MMG::Optimizer().setHMax(hmax)
                     .setHMin(hmin)
-                    .setHausdorff(hausd)
                     .setAngleDetection(false)
                     .optimize(th);
   }