ikarus-project / ikarus

Ikarus is a C++-based library that uses DUNE modules to solve partial differential equations with the finite element method and more.
https://ikarus-project.github.io/
Other
5 stars 3 forks source link

Write a wrapper to fix Dirichlet BCs for each Lagrange nodes #180

Closed tarun-mitruka closed 8 months ago

tarun-mitruka commented 1 year ago

Example:

///// Dune Book Page 314 - for 3D (hard-coded)
template <typename ElementType, typename LocalFE>
void obtainLagrangeNodePositions(const ElementType &ele, const LocalFE &localFE,
                                 std::vector<Dune::FieldVector<double, 3>> &lagrangeNodeCoords) {
  lagrangeNodeCoords.resize(localFE.size());
  std::vector<double> out;
  for (int i = 0; i < 3; i++) {
    auto ithCoord = [&i](const Dune::FieldVector<double, 3> &x) { return x[i]; };
    localFE.localInterpolation().interpolate(ithCoord, out);
    for (std::size_t j = 0; j < out.size(); j++)
      lagrangeNodeCoords[j][i] = out[j];
  }
  for (auto &nCoord : lagrangeNodeCoords)
    nCoord = ele.geometry().global(nCoord);
}

Ikarus::DirichletValues dirichletValues(basisP->flat());

  auto localView = basis.flat().localView();
  for (auto &ele : elements(gridView)) {
    localView.bind(ele);
    auto &fe = localView.tree().child(_0).finiteElement();
    std::vector<Dune::FieldVector<double, 3>> nodalPos;
    obtainLagrangeNodePositions(ele, fe, nodalPos);
    auto tol = 1e-8;
    for (int i = 0; i < fe.size(); i++) {
      if (std::abs(nodalPos[i][0]) < tol) {
        auto fixIndex = localView.index(localView.tree().child(_0).localIndex(i));
        dirichletValues.fixIthDOF(fixIndex);
      }
      if (std::abs(nodalPos[i][1]) < tol) {
        auto fixIndex = localView.index(localView.tree().child(_1).localIndex(i));
        dirichletValues.fixIthDOF(fixIndex);
      }
    }
  }

In dirichletValues.hh,

void fixIthDOF(std::size_t i){dirichletFlags[i] = true;}