path: root/extern/mantaflow/preprocessed/plugin/vortexplugins.cpp

// DO NOT EDIT !
// This file is generated using the MantaFlow preprocessor (prep generate).

/******************************************************************************
 *
 * MantaFlow fluid solver framework
 * Copyright 2011 Tobias Pfaff, Nils Thuerey
 *
 * This program is free software, distributed under the terms of the
 * Apache License, Version 2.0
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Plugins for using vortex sheet meshes
 *
 ******************************************************************************/

#include <iostream>
#include "vortexsheet.h"
#include "vortexpart.h"
#include "shapes.h"
#include "commonkernels.h"
#include "conjugategrad.h"
#include "randomstream.h"
#include "levelset.h"

using namespace std;

namespace Manta {

//! Mark area of mesh inside shape as fixed nodes.
//! Remove all other fixed nodes if 'exclusive' is set

void markAsFixed(Mesh &mesh, const Shape *shape, bool exclusive = true)
{
  for (int i = 0; i < mesh.numNodes(); i++) {
    if (shape->isInside(mesh.nodes(i).pos))
      mesh.nodes(i).flags |= Mesh::NfFixed;
    else if (exclusive)
      mesh.nodes(i).flags &= ~Mesh::NfFixed;
  }
}
static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "markAsFixed", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      Mesh &mesh = *_args.getPtr<Mesh>("mesh", 0, &_lock);
      const Shape *shape = _args.getPtr<Shape>("shape", 1, &_lock);
      bool exclusive = _args.getOpt<bool>("exclusive", 2, true, &_lock);
      _retval = getPyNone();
      markAsFixed(mesh, shape, exclusive);
      _args.check();
    }
    pbFinalizePlugin(parent, "markAsFixed", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("markAsFixed", e.what());
    return 0;
  }
}
static const Pb::Register _RP_markAsFixed("", "markAsFixed", _W_0);
extern "C" {
void PbRegister_markAsFixed()
{
  KEEP_UNUSED(_RP_markAsFixed);
}
}

//! Adapt texture coordinates of mesh inside shape
//! to obtain an effective inflow effect

void texcoordInflow(VortexSheetMesh &mesh, const Shape *shape, const MACGrid &vel)
{
  static Vec3 t0 = Vec3::Zero;

  // get mean velocity
  int cnt = 0;
  Vec3 meanV(0.0);
  FOR_IJK(vel)
  {
    if (shape->isInsideGrid(i, j, k)) {
      cnt++;
      meanV += vel.getCentered(i, j, k);
    }
  }
  meanV /= (Real)cnt;
  t0 -= mesh.getParent()->getDt() * meanV;
  mesh.setReferenceTexOffset(t0);

  // apply mean velocity
  for (int i = 0; i < mesh.numNodes(); i++) {
    if (shape->isInside(mesh.nodes(i).pos)) {
      Vec3 tc = mesh.nodes(i).pos + t0;
      mesh.tex1(i) = tc;
      mesh.tex2(i) = tc;
    }
  }
}
static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "texcoordInflow", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
      const Shape *shape = _args.getPtr<Shape>("shape", 1, &_lock);
      const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 2, &_lock);
      _retval = getPyNone();
      texcoordInflow(mesh, shape, vel);
      _args.check();
    }
    pbFinalizePlugin(parent, "texcoordInflow", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("texcoordInflow", e.what());
    return 0;
  }
}
static const Pb::Register _RP_texcoordInflow("", "texcoordInflow", _W_1);
extern "C" {
void PbRegister_texcoordInflow()
{
  KEEP_UNUSED(_RP_texcoordInflow);
}
}

;

//! Init smoke density values of the mesh surface inside source shape

void meshSmokeInflow(VortexSheetMesh &mesh, const Shape *shape, Real amount)
{
  for (int t = 0; t < mesh.numTris(); t++) {
    if (shape->isInside(mesh.getFaceCenter(t)))
      mesh.sheet(t).smokeAmount = amount;
  }
}
static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "meshSmokeInflow", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
      const Shape *shape = _args.getPtr<Shape>("shape", 1, &_lock);
      Real amount = _args.get<Real>("amount", 2, &_lock);
      _retval = getPyNone();
      meshSmokeInflow(mesh, shape, amount);
      _args.check();
    }
    pbFinalizePlugin(parent, "meshSmokeInflow", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("meshSmokeInflow", e.what());
    return 0;
  }
}
static const Pb::Register _RP_meshSmokeInflow("", "meshSmokeInflow", _W_2);
extern "C" {
void PbRegister_meshSmokeInflow()
{
  KEEP_UNUSED(_RP_meshSmokeInflow);
}
}

struct KnAcceleration : public KernelBase {
  KnAcceleration(MACGrid &a, const MACGrid &v1, const MACGrid &v0, const Real idt)
      : KernelBase(&a, 0), a(a), v1(v1), v0(v0), idt(idt)
  {
    runMessage();
    run();
  }
  inline void op(
      IndexInt idx, MACGrid &a, const MACGrid &v1, const MACGrid &v0, const Real idt) const
  {
    a[idx] = (v1[idx] - v0[idx]) * idt;
  }
  inline MACGrid &getArg0()
  {
    return a;
  }
  typedef MACGrid type0;
  inline const MACGrid &getArg1()
  {
    return v1;
  }
  typedef MACGrid type1;
  inline const MACGrid &getArg2()
  {
    return v0;
  }
  typedef MACGrid type2;
  inline const Real &getArg3()
  {
    return idt;
  }
  typedef Real type3;
  void runMessage()
  {
    debMsg("Executing kernel KnAcceleration ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, a, v1, v0, idt);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MACGrid &a;
  const MACGrid &v1;
  const MACGrid &v0;
  const Real idt;
};

//! Add vorticity to vortex sheets based on buoyancy

void vorticitySource(VortexSheetMesh &mesh,
                     Vec3 gravity,
                     const MACGrid *vel = NULL,
                     const MACGrid *velOld = NULL,
                     Real scale = 0.1,
                     Real maxAmount = 0,
                     Real mult = 1.0)
{
  Real dt = mesh.getParent()->getDt();
  Real dx = mesh.getParent()->getDx();
  MACGrid acceleration(mesh.getParent());
  if (vel)
    KnAcceleration(acceleration, *vel, *velOld, 1.0 / dt);
  const Real A = -1.0;
  Real maxV = 0, meanV = 0;

  for (int t = 0; t < mesh.numTris(); t++) {
    Vec3 fn = mesh.getFaceNormal(t);
    Vec3 source;
    if (vel) {
      Vec3 a = acceleration.getInterpolated(mesh.getFaceCenter(t));
      source = A * cross(fn, a - gravity) * scale;
    }
    else {
      source = A * cross(fn, -gravity) * scale;
    }

    if (mesh.isTriangleFixed(t))
      source = 0;

    mesh.sheet(t).vorticity *= mult;
    mesh.sheet(t).vorticity += dt * source / dx;
    // upper limit
    Real v = norm(mesh.sheet(t).vorticity);
    if (maxAmount > 0 && v > maxAmount)
      mesh.sheet(t).vorticity *= maxAmount / v;

    // stats
    if (v > maxV)
      maxV = v;
    meanV += v;
  }

  cout << "vorticity: max " << maxV << " / mean " << meanV / mesh.numTris() << endl;
}
static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "vorticitySource", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
      Vec3 gravity = _args.get<Vec3>("gravity", 1, &_lock);
      const MACGrid *vel = _args.getPtrOpt<MACGrid>("vel", 2, NULL, &_lock);
      const MACGrid *velOld = _args.getPtrOpt<MACGrid>("velOld", 3, NULL, &_lock);
      Real scale = _args.getOpt<Real>("scale", 4, 0.1, &_lock);
      Real maxAmount = _args.getOpt<Real>("maxAmount", 5, 0, &_lock);
      Real mult = _args.getOpt<Real>("mult", 6, 1.0, &_lock);
      _retval = getPyNone();
      vorticitySource(mesh, gravity, vel, velOld, scale, maxAmount, mult);
      _args.check();
    }
    pbFinalizePlugin(parent, "vorticitySource", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("vorticitySource", e.what());
    return 0;
  }
}
static const Pb::Register _RP_vorticitySource("", "vorticitySource", _W_3);
extern "C" {
void PbRegister_vorticitySource()
{
  KEEP_UNUSED(_RP_vorticitySource);
}
}

void smoothVorticity(VortexSheetMesh &mesh, int iter = 1, Real sigma = 0.2, Real alpha = 0.8)
{
  const Real mult = -0.5 / sigma / sigma;

  // pre-calculate positions and weights
  vector<Vec3> vort(mesh.numTris()), pos(mesh.numTris());
  vector<Real> weights(3 * mesh.numTris());
  vector<int> index(3 * mesh.numTris());
  for (int i = 0; i < mesh.numTris(); i++) {
    pos[i] = mesh.getFaceCenter(i);
    mesh.sheet(i).vorticitySmoothed = mesh.sheet(i).vorticity;
  }
  for (int i = 0; i < mesh.numTris(); i++) {
    for (int c = 0; c < 3; c++) {
      int oc = mesh.corners(i, c).opposite;
      if (oc >= 0) {
        int t = mesh.corners(oc).tri;
        weights[3 * i + c] = exp(normSquare(pos[t] - pos[i]) * mult);
        index[3 * i + c] = t;
      }
      else {
        weights[3 * i + c] = 0;
        index[3 * i + c] = 0;
      }
    }
  }

  for (int it = 0; it < iter; ++it) {
    // first, preload
    for (int i = 0; i < mesh.numTris(); i++)
      vort[i] = mesh.sheet(i).vorticitySmoothed;

    for (int i = 0, idx = 0; i < mesh.numTris(); i++) {
      // loop over adjacent tris
      Real sum = 1.0f;
      Vec3 v = vort[i];
      for (int c = 0; c < 3; c++, idx++) {
        Real w = weights[index[idx]];
        v += w * vort[index[idx]];
        sum += w;
      }
      mesh.sheet(i).vorticitySmoothed = v / sum;
    }
  }
  for (int i = 0; i < mesh.numTris(); i++)
    mesh.sheet(i).vorticitySmoothed *= alpha;
}
static PyObject *_W_4(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "smoothVorticity", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
      int iter = _args.getOpt<int>("iter", 1, 1, &_lock);
      Real sigma = _args.getOpt<Real>("sigma", 2, 0.2, &_lock);
      Real alpha = _args.getOpt<Real>("alpha", 3, 0.8, &_lock);
      _retval = getPyNone();
      smoothVorticity(mesh, iter, sigma, alpha);
      _args.check();
    }
    pbFinalizePlugin(parent, "smoothVorticity", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("smoothVorticity", e.what());
    return 0;
  }
}
static const Pb::Register _RP_smoothVorticity("", "smoothVorticity", _W_4);
extern "C" {
void PbRegister_smoothVorticity()
{
  KEEP_UNUSED(_RP_smoothVorticity);
}
}

//! Seed Vortex Particles inside shape with K41 characteristics
void VPseedK41(VortexParticleSystem &system,
               const Shape *shape,
               Real strength = 0,
               Real sigma0 = 0.2,
               Real sigma1 = 1.0,
               Real probability = 1.0,
               Real N = 3.0)
{
  Grid<Real> temp(system.getParent());
  const Real dt = system.getParent()->getDt();
  static RandomStream rand(3489572);
  Real s0 = pow((Real)sigma0, (Real)(-N + 1.0));
  Real s1 = pow((Real)sigma1, (Real)(-N + 1.0));

  FOR_IJK(temp)
  {
    if (shape->isInsideGrid(i, j, k)) {
      if (rand.getReal() < probability * dt) {
        Real p = rand.getReal();
        Real sigma = pow((1.0 - p) * s0 + p * s1, 1. / (-N + 1.0));
        Vec3 randDir(rand.getReal(), rand.getReal(), rand.getReal());
        Vec3 posUpd(i + rand.getReal(), j + rand.getReal(), k + rand.getReal());
        normalize(randDir);
        Vec3 vorticity = randDir * strength * pow((Real)sigma, (Real)(-10. / 6. + N / 2.0));
        system.add(VortexParticleData(posUpd, vorticity, sigma));
      }
    }
  }
}
static PyObject *_W_5(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "VPseedK41", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      VortexParticleSystem &system = *_args.getPtr<VortexParticleSystem>("system", 0, &_lock);
      const Shape *shape = _args.getPtr<Shape>("shape", 1, &_lock);
      Real strength = _args.getOpt<Real>("strength", 2, 0, &_lock);
      Real sigma0 = _args.getOpt<Real>("sigma0", 3, 0.2, &_lock);
      Real sigma1 = _args.getOpt<Real>("sigma1", 4, 1.0, &_lock);
      Real probability = _args.getOpt<Real>("probability", 5, 1.0, &_lock);
      Real N = _args.getOpt<Real>("N", 6, 3.0, &_lock);
      _retval = getPyNone();
      VPseedK41(system, shape, strength, sigma0, sigma1, probability, N);
      _args.check();
    }
    pbFinalizePlugin(parent, "VPseedK41", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("VPseedK41", e.what());
    return 0;
  }
}
static const Pb::Register _RP_VPseedK41("", "VPseedK41", _W_5);
extern "C" {
void PbRegister_VPseedK41()
{
  KEEP_UNUSED(_RP_VPseedK41);
}
}

//! Vortex-in-cell integration

void VICintegration(VortexSheetMesh &mesh,
                    Real sigma,
                    Grid<Vec3> &vel,
                    const FlagGrid &flags,
                    Grid<Vec3> *vorticity = NULL,
                    Real cgMaxIterFac = 1.5,
                    Real cgAccuracy = 1e-3,
                    Real scale = 0.01,
                    int precondition = 0)
{

  MuTime t0;
  const Real fac = 16.0;  // experimental factor to balance out regularization

  // if no vort grid is given, use a temporary one
  Grid<Vec3> vortTemp(mesh.getParent());
  Grid<Vec3> &vort = (vorticity) ? (*vorticity) : (vortTemp);
  vort.clear();

  // map vorticity to grid using Peskin kernel
  int sgi = ceil(sigma);
  Real pkfac = M_PI / sigma;
  const int numTris = mesh.numTris();
  for (int t = 0; t < numTris; t++) {
    Vec3 pos = mesh.getFaceCenter(t);
    Vec3 v = mesh.sheet(t).vorticity * mesh.getFaceArea(t) * fac;

    // inner kernel
    // first, summate
    Real sum = 0;
    for (int i = -sgi; i < sgi; i++) {
      if (pos.x + i < 0 || (int)pos.x + i >= vort.getSizeX())
        continue;
      for (int j = -sgi; j < sgi; j++) {
        if (pos.y + j < 0 || (int)pos.y + j >= vort.getSizeY())
          continue;
        for (int k = -sgi; k < sgi; k++) {
          if (pos.z + k < 0 || (int)pos.z + k >= vort.getSizeZ())
            continue;
          Vec3i cell(pos.x + i, pos.y + j, pos.z + k);
          if (!flags.isFluid(cell))
            continue;
          Vec3 d = pos -
                   Vec3(i + 0.5 + floor(pos.x), j + 0.5 + floor(pos.y), k + 0.5 + floor(pos.z));
          Real dl = norm(d);
          if (dl > sigma)
            continue;
          // precalc Peskin kernel
          sum += 1.0 + cos(dl * pkfac);
        }
      }
    }
    // then, apply normalized kernel
    Real wnorm = 1.0 / sum;
    for (int i = -sgi; i < sgi; i++) {
      if (pos.x + i < 0 || (int)pos.x + i >= vort.getSizeX())
        continue;
      for (int j = -sgi; j < sgi; j++) {
        if (pos.y + j < 0 || (int)pos.y + j >= vort.getSizeY())
          continue;
        for (int k = -sgi; k < sgi; k++) {
          if (pos.z + k < 0 || (int)pos.z + k >= vort.getSizeZ())
            continue;
          Vec3i cell(pos.x + i, pos.y + j, pos.z + k);
          if (!flags.isFluid(cell))
            continue;
          Vec3 d = pos -
                   Vec3(i + 0.5 + floor(pos.x), j + 0.5 + floor(pos.y), k + 0.5 + floor(pos.z));
          Real dl = norm(d);
          if (dl > sigma)
            continue;
          Real w = (1.0 + cos(dl * pkfac)) * wnorm;
          vort(cell) += v * w;
        }
      }
    }
  }

  // Prepare grids for poisson solve
  Grid<Vec3> vortexCurl(mesh.getParent());
  Grid<Real> rhs(mesh.getParent());
  Grid<Real> solution(mesh.getParent());
  Grid<Real> residual(mesh.getParent());
  Grid<Real> search(mesh.getParent());
  Grid<Real> temp1(mesh.getParent());
  Grid<Real> A0(mesh.getParent());
  Grid<Real> Ai(mesh.getParent());
  Grid<Real> Aj(mesh.getParent());
  Grid<Real> Ak(mesh.getParent());
  Grid<Real> pca0(mesh.getParent());
  Grid<Real> pca1(mesh.getParent());
  Grid<Real> pca2(mesh.getParent());
  Grid<Real> pca3(mesh.getParent());

  MakeLaplaceMatrix(flags, A0, Ai, Aj, Ak);
  CurlOp(vort, vortexCurl);

  // Solve vector poisson equation
  for (int c = 0; c < 3; c++) {
    // construct rhs
    if (vel.getType() & GridBase::TypeMAC)
      GetShiftedComponent(vortexCurl, rhs, c);
    else
      GetComponent(vortexCurl, rhs, c);

    // prepare CG solver
    const int maxIter = (int)(cgMaxIterFac * vel.getSize().max());
    GridCgInterface *gcg = new GridCg<ApplyMatrix>(
        solution, rhs, residual, search, flags, temp1, &A0, &Ai, &Aj, &Ak);
    gcg->setAccuracy(cgAccuracy);
    gcg->setUseL2Norm(true);
    gcg->setICPreconditioner(
        (GridCgInterface::PreconditionType)precondition, &pca0, &pca1, &pca2, &pca3);

    // iterations
    for (int iter = 0; iter < maxIter; iter++) {
      if (!gcg->iterate())
        iter = maxIter;
    }
    debMsg("VICintegration CG iterations:" << gcg->getIterations() << ", res:" << gcg->getSigma(),
           1);
    delete gcg;

    // copy back
    solution *= scale;
    SetComponent(vel, solution, c);
  }
}
static PyObject *_W_6(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "VICintegration", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
      Real sigma = _args.get<Real>("sigma", 1, &_lock);
      Grid<Vec3> &vel = *_args.getPtr<Grid<Vec3>>("vel", 2, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 3, &_lock);
      Grid<Vec3> *vorticity = _args.getPtrOpt<Grid<Vec3>>("vorticity", 4, NULL, &_lock);
      Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 5, 1.5, &_lock);
      Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 6, 1e-3, &_lock);
      Real scale = _args.getOpt<Real>("scale", 7, 0.01, &_lock);
      int precondition = _args.getOpt<int>("precondition", 8, 0, &_lock);
      _retval = getPyNone();
      VICintegration(
          mesh, sigma, vel, flags, vorticity, cgMaxIterFac, cgAccuracy, scale, precondition);
      _args.check();
    }
    pbFinalizePlugin(parent, "VICintegration", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("VICintegration", e.what());
    return 0;
  }
}
static const Pb::Register _RP_VICintegration("", "VICintegration", _W_6);
extern "C" {
void PbRegister_VICintegration()
{
  KEEP_UNUSED(_RP_VICintegration);
}
}

//! Obtain density field from levelset with linear gradient of size sigma over the interface
void densityFromLevelset(const LevelsetGrid &phi,
                         Grid<Real> &density,
                         Real value = 1.0,
                         Real sigma = 1.0)
{
  FOR_IJK(phi)
  {
    // remove boundary
    if (i < 2 || j < 2 || k < 2 || i >= phi.getSizeX() - 2 || j >= phi.getSizeY() - 2 ||
        k >= phi.getSizeZ() - 2)
      density(i, j, k) = 0;
    else if (phi(i, j, k) < -sigma)
      density(i, j, k) = value;
    else if (phi(i, j, k) > sigma)
      density(i, j, k) = 0;
    else
      density(i, j, k) = clamp(
          (Real)(0.5 * value / sigma * (1.0 - phi(i, j, k))), (Real)0.0, value);
  }
}
static PyObject *_W_7(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "densityFromLevelset", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const LevelsetGrid &phi = *_args.getPtr<LevelsetGrid>("phi", 0, &_lock);
      Grid<Real> &density = *_args.getPtr<Grid<Real>>("density", 1, &_lock);
      Real value = _args.getOpt<Real>("value", 2, 1.0, &_lock);
      Real sigma = _args.getOpt<Real>("sigma", 3, 1.0, &_lock);
      _retval = getPyNone();
      densityFromLevelset(phi, density, value, sigma);
      _args.check();
    }
    pbFinalizePlugin(parent, "densityFromLevelset", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("densityFromLevelset", e.what());
    return 0;
  }
}
static const Pb::Register _RP_densityFromLevelset("", "densityFromLevelset", _W_7);
extern "C" {
void PbRegister_densityFromLevelset()
{
  KEEP_UNUSED(_RP_densityFromLevelset);
}
}

}  // namespace Manta