proxpde/main.cpp at master · capitalaslash/proxpde · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
#include "def.hpp"

#include "assembly.hpp"
#include "bc.hpp"
#include "builder.hpp"
#include "fe.hpp"
#include "fespace.hpp"
#include "iomanager.hpp"
#include "mesh.hpp"

enum class SolverType : uint8_t
{
  CHOLESKY,
  BICGSTAB,
  SPARSELU
};

int main()
{
  using namespace proxpde;

  using Elem_T = Quad;
  using Mesh_T = Mesh<Elem_T>;
  using FESpace_T = FESpace<Mesh_T, LagrangeFE<Elem_T, 1>::RefFE_T, GaussQR<Elem_T, 4>>;

  auto const solverType = SolverType::SPARSELU;

  auto const rhs = [](Vec3 const & p) { return M_PI * std::sin(M_PI * p(0)); };

  auto const exactSol = [](Vec3 const & p)
  { return std::sin(M_PI * p(0)) / M_PI + p(0); };

  // auto const rhs = [] (Vec3 const& p) { return p(0); };
  // auto const exact_sol = [] (Vec3 const& p) { return 0.5*p(0) -  p(0)*p(0)*p(0)/6.;
  // };

  // auto const rhs = [] (Vec3 const&) { return 8.; };
  // auto const exact_sol = [] (Vec3 const& p) { return 4.*p(0)*(2.-p(0)); };
  // auto const exact_sol = [] (Vec3 const& p) { return 4.*p(0)*(1.-p(0)); };

  uint const numElemsX = 20;
  uint const numElemsY = 1;

  Vec3 const origin{0., 0., 0.};
  Vec3 const length{1., 0.02, 0.};

  std::unique_ptr<Mesh_T> mesh{new Mesh_T};

  buildHyperCube(*mesh, origin, length, {numElemsX, numElemsY, 0});
  // std::cout << *mesh << std::endl;

  FESpace_T feSpace{*mesh};

  auto bc = BCEss{feSpace, side::LEFT};
  bc << [](Vec3 const &) { return 0.; };

  AssemblyStiffness stiffness{1.0, feSpace};
  AssemblyRhsAnalytic f{rhs, feSpace};

  Builder builder{feSpace.dof.size};
  builder.buildLhs(std::tuple{stiffness}, std::vector{bc});
  builder.buildRhs(std::tuple{f}, std::vector{bc});
  builder.closeMatrix();

  // std::cout << "builder.A:\n" << builder.A << std::endl;
  // std::cout << "b:\n" << builder.b << std::endl;

  // std::ofstream fout("output.m");
  // for( int k=0; k<builder.A.outerSize(); k++)
  // {
  //   for (Mat::InnerIterator it(builder.A,k); it; ++it)
  //   {
  //     std::cout << it.row() << " " << it.col() << " " << it.value() << " " <<
  //     it.index() << std::endl; fout << it.row()+1 << " " << it.col()+1 << " " <<
  //     it.value() << std::endl;
  //   }
  //   std::cout << "-----" << std::endl;
  // }
  // std::cout << "=====" << std::endl;
  // fout.close();

  Var sol{"u"};
  switch (solverType)
  {
  case SolverType::CHOLESKY:
  {
    Eigen::SimplicialCholesky<SparseMatrix<StorageType::ClmMajor>> solver(builder.A);
    sol.data = solver.solve(builder.b);
    break;
  }
  case SolverType::BICGSTAB:
  {
    Eigen::BiCGSTAB<SparseMatrix<StorageType::RowMajor>> solver(builder.A);
    sol.data = solver.solve(builder.b);
    break;
  }
  case SolverType::SPARSELU:
  {
    Eigen::SparseLU<SparseMatrix<StorageType::ClmMajor>, Eigen::COLAMDOrdering<int>>
        solver;
    // Compute the ordering permutation vector from the structural pattern of A
    solver.analyzePattern(builder.A);
    // Compute the numerical factorization
    solver.factorize(builder.A);

    sol.data = solver.solve(builder.b);
    break;
  }
  default:
    std::abort();
  }
  // std::cout<< "sol:\n" << sol << std::endl;

  Var exact{"exact"};
  interpolateAnalyticFunction(exactSol, feSpace, exact.data);

  Var error{"error"};
  error.data = sol.data - exact.data;
  fmt::print("error: {:e}\n", error.data.norm());

  IOManager io{feSpace, "output/sol"};
  io.print({sol, exact, error});

  return 0;
}