Writing zoneID with isoSurface

I'm trying to save an alpha based isoSurface from a case running with overInterDyMFoam.

Everything works fine but if I ask for zoneID as one of the fields to be written with the surface along with velocity and turbulent kinetic energy, only U and k files are produced.

I've tried to put the "dummy" keyword in the "fields" list and I get a warning for it only, so it looks like zoneID is found indeed as registered object.

## Reattaching the author to the issue ticket: @Ricky-11 ##

Probably that we normally sample volScalarField, volVectorField etc, but the zoneID is a labelField (if I remember correctly).

I see. Would it be fairly easy for you to add the labelField to sampling or could you point me in the right direction to try on my own?

By Riccardo Rossi on 2019-07-24T07:19:15 (imported from GitLab project)

## Reattaching the author to the issue ticket: @Ricky-11 ##

It's "medium easy" but I can't look at it for a while (I'm at the OpenFOAM Workshop this week). The sampledSurfaces.C file around line 574 would need an extra template call. But you'd then need to make another implementation of performAction in the Templates file (similar to what is at line 80-170, but dimensionless etc). The field count will still not be quite right, but only affects legacy vtk output.

Another approach, brute force but fast. Have a simple coded function object to create a zoneIds volScalarField from the zoneID labelField by copying and give it a zero gradient bc. Sample the newly defined field. It may not be quite as sleek, but should be quick to implement and get you on the way.

Sounds good. I'll try to take a look to both and I'll keep you posted. Thanks so much for the help.

By Riccardo Rossi on 2019-07-24T08:50:30 (imported from GitLab project)

So I've been successful in creating the new volScalarField with the coded function object but the isoSurface won't find it listed as a registered object for the "fields" list.

I've read on the forums that this might be due to the fact that the new field is visible only within the coded snippet and then destroyed afterwards.

By Riccardo Rossi on 2019-07-24T16:10:34 (imported from GitLab project)

If you create a volScalarField as a pointer (eg, with new) you would use its 'store()' method to flag it as being owned by its registry, which should make it persistant. This is a fairly common approach (sorry no code access at the moment) - see if a vol field exists, if not create one and register+store.

label feature request

I have to correct myself.

The coded function object supposed to create the volScalarField for zoneID wasn't working, apparently because was not at the beginning of the functions list in controlDict (is this mandatory???).

Once properly located, it works but after quite an extensive research in the source code and online I'm still having a hard time to access the zoneID labelField.

Last trial was:

codeExecute #{ volScalarField myZoneID ( IOobject ( "myZoneID", mesh().time().timeName(), mesh(), IOobject::NO_READ, IOobject::AUTO_WRITE ), mesh().lookupObject("zoneID") ); #};

By Riccardo Rossi on 2019-07-25T08:29:30 (imported from GitLab project)

You want a modified logic from src/overset/cellCellStencil/cellCellStencil/cellCellStencil.C (around line 100).

In your case, however, you want to look for zoneID with something like this:

const labelIOList *zoneIDPtr = mesh.cfindObject<labelIOList>("zoneID");

volScalarField* volZonePtr = mesh.getObjectPtr<volScalarField("volZoneID");

if (zoneIDPtr)
{
    if (!volZonePtr)
    {
        volZonePtr = new volScalarField
        (
            IOobject
            (
                "volZoneID",
                mesh.time().findInstance(mesh.dbDir(), "zoneID"),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                true  // Register
             ),
             mesh,
             dimensionedScalar(dimless, Zero),
             zeroGradientFvPatchField<scalar>::typeName
         );

        volZonePtr->store();
    }

    std::copy(zoneIDPtr->begin(), zoneIDPtr->end(), volZonePtr->begin());
}

This is fully untested, but gives the general idea. You may wish to play with the read/write options, depending on your requirements. Could also use a dimensionedSet instead dimensionedScalar in the constructor.

Thanks Mark. I'll give a try right away.

By Riccardo Rossi on 2019-07-25T08:56:22 (imported from GitLab project)

Just tried the code you suggested and got a bunch of error at runtime when the coded object is compiled.

The first one is:

/gpfs/work/rfd_prod1/boardsports/Sequoia/tmp/EnzoOversetNewRunWithFreeSurfaceClipping/system/controlDict.functions.myZoneID:88:32: error: ‘((Foam::myZoneIDFunctionObject*)this)->Foam::myZoneIDFunctionObject::mesh’ does not have class type

If you have the chance to take a look at this, it might be helpful for issue #1387 to compute and write wall distance along with the isosurface.

By Riccardo Rossi on 2019-08-03T21:08:20 (imported from GitLab project)

Also tried to work on this but no luck. Any additional guidance @mark?

By Riccardo Rossi on 2019-08-20T21:59:57 (imported from GitLab project)

Just adding an example now, will reference this issue and push in a bit (building/testing some other things first).

mentioned in commit d043850bbcec2463f4aa2b5622b386a703ed7696

Didn't get any notification by mail of your commit, so I only found out today by checking the devel site.

The example works silky smooth with the angledDuct/implicit as well as with the twoSimpleRotors tutorial.

Unfortunately, the code won't compile with the v1706, which is the only working version for our problem at the moment (see #1183 (closed)).

I'll take the chance to go back to the other issue and see if we can finally solve it with the help of @Mattijs

By Riccardo Rossi on 2019-08-26T16:02:06 (imported from GitLab project)

Hi @Ricky-11 ,

Didn't figure on seeing 1706 again, but dug out a version. As suspected, there are too many changes needed to demodernize the code. Try the following:

Instead of the shorter names,

libs    ("utilityFunctionObjects");

you'll need the longer ones (for 1812 and earlier),

libs    ("libutilityFunctionObjects.so");
libs    ("libsampling.so");

You will also need the older (IMO more confusing lookup name). Instead of

volScalarField* volZonePtr = mesh.getObjectPtr<volScalarField>(fieldName);

try with

volScalarField* volZonePtr = mesh.lookupObjectRefPtr<volScalarField>(fieldName);

Finally, you will need to provide the interpolationScheme for older versions:

interpolationScheme cellPoint;

Thanks for looking into the v1706 version @mark.

I had figure out right away libs calls had changed, but then I got stuck with the pointer, so I moved on and fixed the case setup in the v1906.

The use of the zoneID as a marker to cleanup the isoSurface within the overInterDyMFoam run didn't go as expected, though.

Attached are two pictures showing a comparison between both the isoSurface and cutPlane clipped with zone id in Paraview starting from reading the case and from the vtk surfaces carrying the zoneID.

As you can see, the isoSurface created in paraview shows no interpolation of the zoneid scalar between the background and the moving mesh and once clipped only the trianlges lying in the background mesh are kept and no holes are present.

On the other hand, the zone id read with the vtk surface written at runtime is forcefully interpolated (of complains about the need to cross cells if interpolation is deactivated) and once clipped the triangles with the intermediate zone id values remains plus holes in the background mesh are found.

Even with cutplanes, however, where interpolation is set to off and the zone id separation is crisp, once the plane is clipped holes in the background mesh are also found.

Any thought?

By Riccardo Rossi on 2019-08-26T23:12:17 (imported from GitLab project)

Any chance to get any hint/direction on this @mark?

By Riccardo Rossi on 2019-10-10T17:08:35 (imported from GitLab project)

changed the description

Closing ticket - please re-open if still relevant

closed