Merge pull request #99 from cbritopacheco/develop

Complex Numbers

examples/BoundaryOptimization/CathodeAnode.cpp

@@ -172,8 +172,6 @@ int main(int, char**)
     // dOmega.save("out/Anode." + std::to_string(i) + ".mesh");
     // distAnode.save("out/Anode." + std::to_string(i) + ".gf");
 
-    Solver::CG cg;
-
     // Parameters
     RealFunction uIn = 1.0;
     RealFunction gamma = 1.0;
@@ -201,7 +199,7 @@ int main(int, char**)
           - BoundaryIntegral(heAnode * uIn, v)
           ;
     Alert::Info() << "Solving state equation..." << Alert::Raise;
-    state.solve(cg);
+    Solver::CG(state).solve();
     u.getSolution().save("U.gf");
     mesh.save("Omega.mfem.mesh");
 
@@ -212,7 +210,7 @@ int main(int, char**)
               + Integral(gtr, gte)
               - Integral(0.5 * Grad(u.getSolution()) / mesh.getVolume(), gte)
               ;
-    potential.solve(cg);
+    Solver::CG(potential).solve();
 
     Alert::Info() << "Solving adjoint equation..." << Alert::Raise;
     TrialFunction p(sfes);
@@ -221,7 +219,7 @@ int main(int, char**)
     adjoint = Integral(gamma * Grad(p), Grad(q))
             + BoundaryIntegral((heAnode + heCathode) * p, q)
             - Integral(gtr.getSolution(), Grad(q));
-    adjoint.solve(cg);
+    Solver::CG(adjoint).solve();
     p.getSolution().save("P.gf");
 
     GridFunction sC(sfes);
@@ -273,7 +271,7 @@ int main(int, char**)
             + lambdaCathode * FaceIntegral(
                 Div(conormalCathode).traceOf(Cathode) * conormalCathode, w).over(dCathode)
             + tgv * Integral(theta, w).over(Inlet, FixedAnode, FixedCathode);
-    hilbert.solve(cg);
+    Solver::CG(hilbert).solve();
 
     GridFunction norm(dsfes);
     norm = Frobenius(theta.getSolution());
@@ -306,7 +304,7 @@ int main(int, char**)
              + Integral((u.getSolution()) * p.getSolution(), t)
              + tgv * Integral(s, t).over(Inlet, Cathode, Anode, FixedAnode, FixedCathode);
       }
-      topo.solve(cg);
+      Solver::CG(topo).solve();
 
       s.getSolution().save("Topo.sol", IO::FileFormat::MEDIT);
       dsfes.getMesh().save("Topo.mesh", IO::FileFormat::MEDIT);

examples/BoundaryOptimization/ClampingLocator.cpp

@@ -269,8 +269,6 @@ int main(int, char**)
       MMG::Distancer(dsfes).setInteriorDomain(Locator)
                            .distance(dOmega);
 
-    Solver::CG cg;
-
     // Parameters
     const Real lambda = 0.5769, mu = 0.3846;
 
@@ -311,7 +309,7 @@ int main(int, char**)
           //+ BoundaryIntegral(tgv * u.z(), v.z()).over(7);;
     state.assemble();
     Alert::Info() << "Solving state equation..." << Alert::Raise;
-    state.solve(cg);
+    Solver::CG(state).solve();
 
     u.getSolution().save("u.gf");
     mesh.save("u.mesh");
@@ -323,8 +321,7 @@ int main(int, char**)
     adjoint = LinearElasticityIntegral(p, q)(lambda, mu)
             + Integral(u.getSolution() / mesh.getVolume(), q)
             + FaceIntegral(heClamp * p, q).over(Clamp);
-    adjoint.solve(cg);
-
+    Solver::CG(adjoint).solve();
 
     p.getSolution().save("p.gf");
     mesh.save("p.mesh");
@@ -374,7 +371,7 @@ int main(int, char**)
              - Integral(Dot(f, p.getSolution()), t)
              + tgv * Integral(s, t).over(Clamp, Locator, GammaT);
       }
-      topo.solve(cg);
+      Solver::CG(topo).solve();
 
       s.getSolution().save("Topo.gf");
       dsfes.getMesh().save("Topo.mesh");
@@ -430,7 +427,7 @@ int main(int, char**)
               + ellClamp * FaceIntegral(conormalClamp, w).over(dClamp)
               + ellLocator * FaceIntegral(conormalLocator, w).over(dLocator)
               + tgv * Integral(theta, w).over(GammaT);
-      hilbert.solve(cg);
+      Solver::CG(hilbert).solve();
 
       GridFunction norm(dsfes);
       norm = Frobenius(theta.getSolution());

examples/BoundaryOptimization/CoolingMaterial.cpp

@@ -143,8 +143,6 @@ int main(int, char**)
                                      .distance(dOmega);
 
     // Parameters
-    Solver::CG cg;
-
     RealFunction f = 1;
     RealFunction g = -1.0;
 
@@ -168,7 +166,7 @@ int main(int, char**)
     state = Integral(Grad(u), Grad(v))
           + FaceIntegral(he * u, v).over({Gamma, GammaD})
           - Integral(f, v);
-    state.solve(cg);
+    Solver::CG(state).solve();
 
     auto dj = -1.0 / Omega.getVolume();
     Alert::Info() << "Solving adjoint equation..." << Alert::Raise;
@@ -178,7 +176,7 @@ int main(int, char**)
     adjoint = Integral(Grad(p), Grad(q))
             + BoundaryIntegral(he * p, q).over({Gamma, GammaD})
             - Integral(dj * q);
-    adjoint.solve(cg);
+    Solver::CG(adjoint).solve();
 
     Alert::Info() << "Computing objective..." << Alert::Raise;
     const Real objective = J(u.getSolution());
@@ -200,7 +198,7 @@ int main(int, char**)
             + Integral(theta, w)
             + Integral(tgv * g, w).over(GammaN)
             - FaceIntegral(hadamard, w).over(SigmaD);
-    hilbert.solve(cg);
+    Solver::CG(hilbert).solve();
     GridFunction grad(dvfes);
     grad = theta.getSolution() * conormal;
     grad *= -1.0;

examples/BoundaryOptimization/WaterTankSupport.cpp

@@ -274,8 +274,6 @@ int main(int, char**)
     auto dist = MMG::Distancer(dsfes).setInteriorDomain(Support)
                                      .distance(dOmega);
 
-    Solver::CG cg;
-
     // Parameters
     const Real lambda = 0.5769, mu = 0.3846;
 
@@ -303,7 +301,7 @@ int main(int, char**)
           - Integral(g, v)
           + FaceIntegral(he * u, v).over(Support);
           ;
-    state.solve(cg);
+    Solver::CG(state).solve();
 
     u.getSolution().save("u.gf");
     mesh.save("u.mesh");
@@ -315,7 +313,7 @@ int main(int, char**)
     adjoint = LinearElasticityIntegral(p, q)(lambda, mu)
             + Integral(u.getSolution() / mesh.getVolume(), q)
             + FaceIntegral(he * p, q).over(Support);
-    adjoint.solve(cg);
+    Solver::CG(adjoint).solve();
 
     p.getSolution().save("p.gf");
     mesh.save("p.mesh");
@@ -358,7 +356,7 @@ int main(int, char**)
            + Integral(s, t)
            + Integral((u.getSolution().T() * aniso * p.getSolution()).coeff(0, 0), t)
            + tgv * Integral(s, t).over(Support, Unsupported);
-      topo.solve(cg);
+      Solver::CG(topo).solve();
 
       s.getSolution().save("Topo.gf");
       dsfes.getMesh().save("Topo.mesh");
@@ -413,7 +411,7 @@ int main(int, char**)
               // + ellP * FaceIntegral(Div(conormal).traceOf(Support) * conormal, w).over(dSupport)
               // + tgv * Integral(theta, w).over(Unsupported)
               ;
-      hilbert.solve(cg);
+      Solver::CG(hilbert).solve();
 
       GridFunction norm(dsfes);
       norm = Frobenius(theta.getSolution());

examples/BoundaryOptimization/WaveConcealing.cpp

@@ -160,8 +160,6 @@ int main(int, char**)
     dist.save("Dist.gf");
     dist.getFiniteElementSpace().getMesh().save("Dist.mesh");
 
-    Solver::CG cg;
-
     // Parameters
     VectorFunction xi = { 0, 0, 1 };
     RealFunction phi =
@@ -197,7 +195,7 @@ int main(int, char**)
                 - Integral(waveNumber * waveNumber * gamma * u0, v0)
                 + DirichletBC(u0, phi).on(PlaneWave, Dirichlet);
                 //+ DirichletBC(u0, Zero()).on(Dirichlet);
-    unperturbed.solve(cg);
+    Solver::CG(unperturbed).solve();
 
     u0.getSolution().save("U0.gf");
     u0.getSolution().getFiniteElementSpace().getMesh().save("U0.mesh");
@@ -211,7 +209,7 @@ int main(int, char**)
           + FaceIntegral(he * u, v).over({ Absorbing, Reflecting, Window })
           + DirichletBC(u, phi).on(PlaneWave, Dirichlet);
           //+ DirichletBC(u, Zero()).on(Dirichlet);
-    state.solve(cg);
+    Solver::CG(state).solve();
 
     u.getSolution().save("State.gf");
     u.getSolution().getFiniteElementSpace().getMesh().save("State.mesh");
@@ -232,7 +230,7 @@ int main(int, char**)
             + FaceIntegral(he * p, q).over({ Absorbing, Reflecting, Window })
             + DirichletBC(u, Zero()).on({ Dirichlet, PlaneWave });
             ;
-    adjoint.solve(cg);
+    Solver::CG(adjoint).solve();
 
     p.getSolution().save("Adjoint.gf");
     p.getSolution().getFiniteElementSpace().getMesh().save("Adjoint.mesh");
@@ -272,7 +270,7 @@ int main(int, char**)
            + Integral(s, t)
            + Integral(u.getSolution() * p.getSolution(), t)
            + tgv * Integral(s, t).over(Absorbing, Window, Dirichlet, PlaneWave);
-      topo.solve(cg);
+      Solver::CG(topo).solve();
 
       s.getSolution().save("Topo.gf");
       dsfes.getMesh().save("Topo.mesh");
@@ -335,7 +333,7 @@ int main(int, char**)
               + 0.5 * 2 * bA * constraint * FaceIntegral(conormal, w).over(dAbsorbing)
               + tgv * Integral(theta, w).over(Dirichlet, PlaneWave, Window)
               ;
-      hilbert.solve(cg);
+      Solver::CG(hilbert).solve();
 
       GridFunction norm(dsfes);
       norm = Frobenius(theta.getSolution());

examples/DensityOptimization/TemperatureMinimization.cpp

@@ -52,15 +52,13 @@ int main(int, char**)
     // Poisson problem
     RealFunction f(1.0);
 
-    Solver::SparseLU solver;
-
     TrialFunction u(vh);
     TestFunction  v(vh);
     Problem poisson(u, v);
     poisson = Integral((gmin + (gmax - gmin) * Pow(gamma, 3)) * Grad(u), Grad(v))
             - Integral(f * v)
             + DirichletBC(u, RealFunction(0.0)).on(GammaD);
-    poisson.solve(solver);
+    Solver::SparseLU(poisson).solve();
 
     // Adjoint problem
     TrialFunction p(vh);
@@ -69,7 +67,7 @@ int main(int, char**)
     adjoint = Integral((gmin + (gmax - gmin) * Pow(gamma, 3)) * Grad(p), Grad(q))
             + Integral(RealFunction(1.0 / vol), q)
             + DirichletBC(p, RealFunction(0.0)).on(GammaD);
-    adjoint.solve(solver);
+    Solver::SparseLU(adjoint).solve();
 
     // Hilbert extension-regularization
     TrialFunction g(vh);
@@ -81,7 +79,7 @@ int main(int, char**)
                 ell + 3 * (gmax - gmin) * Pow(gamma, 2) *
                 Dot(Grad(u.getSolution()), Grad(p.getSolution())), w)
             + DirichletBC(g, RealFunction(0.0)).on(GammaD);
-    hilbert.solve(solver);
+    Solver::SparseLU(hilbert).solve();
 
     GridFunction step(ph);
     step = mu * g.getSolution();

examples/IntegralEquations/ElasticDistribution.cpp

@@ -70,8 +70,7 @@ int main(int, char**)
   // file.close();
 
   std::cout << "resolution\n";
-  Solver::CG<Math::Matrix<Real>, Math::Vector<Real>> cg;
-  eq.solve(cg);
+  Solver::CG(eq).solve();
 
   u.getSolution().save("u.gf");
   mesh.save("u.mesh");

examples/IntegralEquations/NewtonianPotential.cpp

@@ -54,8 +54,7 @@ int main(int, char**)
   file.close();
 
   std::cout << "resolution\n";
-  Solver::LDLT solver;
-  eq.solve(solver);
+  Solver::LDLT(eq).solve();
 
   u.getSolution().save("u.gf");
   mesh.save("u.mesh");

examples/Misc/RS2023/CellFunction.cpp

@@ -83,16 +83,13 @@ int main(int, char**)
            + Integral(dygamma, v)
            + PeriodicBC(psi2, dofs);
 
-  // Solve the problem
-  Solver::SparseLU solver;
-
   Alert::Info() << "Solving for first component..." << Alert::Raise;
-  poisson1.solve(solver);
+  Solver::SparseLU(poisson1).solve();
   psi1.getSolution().save("CellFunction1.gf");
   Alert::Success() << "Done! Saved to CellFunction1.gf" << Alert::Raise;
 
   Alert::Info() << "Solving for second component..." << Alert::Raise;
-  poisson2.solve(solver);
+  Solver::SparseLU(poisson2).solve();
   psi2.getSolution().save("CellFunction2.gf");
   Alert::Success() << "Done! Saved to CellFunction2.gf" << Alert::Raise;
 

examples/Misc/RS2023/Conductivity/Perturbed.cpp

@@ -12,6 +12,7 @@
 
 using namespace Rodin;
 using namespace Rodin::Math;
+using namespace Rodin::Solver;
 using namespace Rodin::Geometry;
 using namespace Rodin::Variational;
 
@@ -36,8 +37,6 @@ static constexpr Real gamma_ek = 1e+12;
 static constexpr Real R0 = 0.2; // Radius of B_R(x_0)
 static constexpr Real R1 = R0 + 10 * hmax; // Radius of B_R(x_0)
 
-static Solver::SparseLU solver;
-
 int main(int, char**)
 {
   Alert::Info() << "Epsilon: " << epsilon
@@ -98,7 +97,7 @@ int main(int, char**)
 
   // Solve the background problem
   Alert::Info() << "Solving background equation." << Alert::Raise;
-  poisson.solve(solver);
+  SparseLU(poisson).solve();
   const auto u0 = u.getSolution();
 
   u0.save("Background.gf");
@@ -111,7 +110,7 @@ int main(int, char**)
   // Solve the perturbed problem
   Alert::Info() << "Solving perturbed equation." << Alert::Raise;
 
-  perturbed.solve(solver);
+  SparseLU(perturbed).solve();
   const auto u_e = u.getSolution();
   GridFunction g_e(gh);
   u_e.save("Perturbed.gf");

examples/Misc/RS2023/Conductivity/Simulation.cpp

@@ -27,8 +27,6 @@ static const Math::Vector<Real> x0{{0.5, 0.5}}; // Center of domain
 static constexpr Real pi = Math::Constants::pi();
 static constexpr Real gamma_ek = 1.93559;
 
-static Solver::SparseLU solver;
-
 struct Data
 {
   Real m;
@@ -140,7 +138,7 @@ int main(int, char**)
 
       // Solve the background problem
       Alert::Info() << "Solving background equation." << Alert::Raise;
-      poisson.assemble().solve(solver);
+      Solver::SparseLU(poisson).solve();
       const auto u0 = u.getSolution();
 
       // u0.save("Background.gf");
@@ -153,7 +151,7 @@ int main(int, char**)
       // Solve the perturbed problem
       Alert::Info() << "Solving perturbed equation." << Alert::Raise;
 
-      perturbed.assemble().solve(solver);
+      Solver::SparseLU(perturbed).solve();
       const auto u_e = u.getSolution();
       // GridFunction g_e(gh);
       // // u_e.save("Perturbed.gf");

examples/Misc/RS2023/Helmholtz/Animation.cpp

@@ -82,8 +82,6 @@ void run(int id, const std::vector<Data>& grid)
   mesh.load("Q1.medit.mesh", IO::FileFormat::MEDIT);
   mesh.save("out/Q.mesh");
 
-  Solver::SparseLU solver;
-
   Alert::Info() << "Grid size: " << grid.size() << Alert::Raise;
 
   size_t i = 0;
@@ -155,12 +153,12 @@ void run(int id, const std::vector<Data>& grid)
 
     // Solve the background problem
     Alert::Info() << "Solving background equation." << Alert::Raise;
-    helmholtz.solve(solver);
+    Solver::SparseLU(helmholtz).solve();
     const auto u0 = std::move(u.getSolution());
     Alert::Success() << "Done." << Alert::Raise;
 
     Alert::Info() << "Solving perturbed equation." << Alert::Raise;
-    perturbed.solve(solver);
+    Solver::SparseLU(perturbed).solve();
     const auto ue = std::move(u.getSolution());
     Alert::Success() << "Done." << Alert::Raise;