test.cpp

#include "Time.H"
#include "argList.H"
#include "fvMesh.H"
#include "volFields.H"
#include "pointFields.H"
#include "dimensionedScalar.H"
#include "constants.H"
#include "interpolatePointToCell.H"

using namespace Foam;
using Foam::constant::mathematical::pi;

int main(int argc, char **argv) {
  argList args(argc, argv);
  if (not args.checkRootCase()) FatalError.exit();
  Time run_time(Time::controlDictName, args);
  fvMesh mesh(IOobject(fvMesh::defaultRegion, run_time.timeName(), run_time,
                       IOobject::MUST_READ));
  const pointMesh& pmesh = pointMesh::New(mesh);

  pointScalarField pf(IOobject("pf", run_time.timeName(), mesh,
                               IOobject::NO_READ, IOobject::AUTO_WRITE),
                      pmesh, dimensionedScalar("0", dimless, 0.0));

  volScalarField vf(IOobject("vf", run_time.timeName(), mesh,
                               IOobject::NO_READ, IOobject::AUTO_WRITE),
                      mesh, dimensionedScalar("0", dimless, 0.0));

  const auto& bb = mesh.bounds();
  const auto& c = bb.midpoint();
  const scalar cx = c.x();
  const scalar cy = c.y();
  forAll(pf, i) {
    const auto& p = mesh.points()[i];
    const scalar r = Foam::sqrt(sqr(p.x() - cx) + sqr(p.y() - cy));
    pf[i] = 2*Foam::j1(pi*r)/(pi*r + SMALL);
  }

  forAll(vf, i) {
    vf[i] = interpolatePointToCell(pf, i);
  }

  pf.write();
  vf.write();

  return 0;
}