[ganeshv]

hi !

I need to perform a mapping from a fine mesh to a coarse mesh every few time steps. The coarse mesh can be a sub-domain of the full fine mesh. This needs to run in parallel. My search for such a utility took me here

http://www.cfd-online.com/Forums/ope...time-step.html

But it's not very conclusive. I thought I could modify the source code of mapFields to do something in this regard. The idea is to read in every processor sub-mesh of the coarse mesh on all the processors. Initiate the desired fields to zero. Map the fields from the fine mesh to the coarse mesh and sum them all up on one processor and write it out. This would essentially become some sort of parallel mapFields utility. I can't get all processors to read in a processor mesh like so

        Time runTimeCoarseMesh
          (
           Time::controlDictName,
           word("."),
           word(".")/fileName(word("processor")+name(procI))
           );

        fvMesh coarseMesh
          (
           IOobject
           (
            "coarseSolution",
            runTimeCoarseMesh.timeName(),
            runTimeCoarseMesh,
            IOobject::MUST_READ
            )
           );

[carowjp]

Hello,

Instead of needing to map from a course to a fine mesh I am mapping from a rotated to non-rotated mesh:

meshToMesh mtmCon(mesh, meshR);
mtmCon.interpolate(alphaR, alpha, meshToMesh::MAP); // map consistent

quaternion Rotation(axis,theta);
tensor R = Rotation.R();
meshPoints = transform(R, meshPoints);
meshR.movePoints(meshPoints);              // rotate the mesh

meshToMesh mtmInc(meshR, mesh, patchMap, cuttingPatches);
mtmInc.interpolate(alpha, alphaR, meshToMesh::INTERPOLATE); // interpolate inconsistent

The cutting patches are read from a mapFieldsDict.

In any case using meshToMesh interpolation in parallel gives incorrect results as shown in the attached image.

The image is a composite of three frames for a case where a scalar field is mapped. The left image is the original, the middle image is rotated and mapped with meshToMesh interpolation in serial (correct), and the right image mapped with meshToMesh interpolation in parallel with the mesh decomposed for 4 processors (2 x 2).

In the parallel case the mapping only takes place on processor0, so the field is 'clipped' by the processor patches.

Can anyone please shed any light on how to address this issue?

thank you,

James

[carowjp]

Sorry for the bump but I am really stumped!

Can anyone please help?

regards,

James

[carowjp]

Hello again,

The meshToMesh interpolation seems to only map on processor0 when run in parallel. I have studied the mapFields utility and attempted to implement the parallel specific portions of code into my application.

With mapFields, the decomposed source and target meshes are read from file:

    
if (parallelSource && parallelTarget)
{
        //after reading source and target "decomposeParDict" files & determining number of subdomains
        
        List<boundBox> bbsTarget(nProcsTarget);
        List<bool> bbsTargetSet(nProcsTarget, false);

       // loop over processors for source mesh
        for (int procISource=0; procISource<nProcsSource; procISource++)
        {
            Info<< nl << "Source processor " << procISource << endl;

            Time runTimeSource
            (
                Time::controlDictName,
                rootDirSource,
                caseDirSource/fileName(word("processor") + name(procISource))
            );

            #include "setTimeIndex.H"

            fvMesh meshSource
            (
                IOobject
                (
                    sourceRegion,
                    runTimeSource.timeName(),
                    runTimeSource
                )
            );

            Info<< "mesh size: " << meshSource.nCells() << endl;

            boundBox bbSource(meshSource.bounds());

            // loop over processors for target mesh
            for (int procITarget=0; procITarget<nProcsTarget; procITarget++)
            {
                if
                (
                    !bbsTargetSet[procITarget]
                  || (
                      bbsTargetSet[procITarget]
                   && bbsTarget[procITarget].overlaps(bbSource)
                     )
                )
                {
                    Info<< nl << "Target processor " << procITarget << endl;

                    Time runTimeTarget
                    (
                        Time::controlDictName,
                        rootDirTarget,
                        caseDirTarget/fileName(word("processor")
                      + name(procITarget))
                    );

                    fvMesh meshTarget
                    (
                        IOobject
                        (
                            targetRegion,
                            runTimeTarget.timeName(),
                            runTimeTarget
                        )
                    );

                    Info<< "mesh size: " << meshTarget.nCells() << endl;

                    bbsTarget[procITarget] = meshTarget.bounds();
                    bbsTargetSet[procITarget] = true;

                    if (bbsTarget[procITarget].overlaps(bbSource))
                    {
                        if (consistent)
                        {
                            mapConsistentSubMesh
                            (
                                meshSource,
                                meshTarget,
                                mapOrder
                            );
                        }
                        else
                        {
                            mapSubMesh
                            (
                                meshSource,
                                meshTarget,
                                patchMap,
                                addProcessorPatches(meshTarget, cuttingPatches),
                                mapOrder
                            );
                        }
                    }
                }
            }
        }
    }
}

Can anyone please explain how to do this at runtime (without reading from file) by addressing the source and target meshes on each processor?

This was discussed in the thread: Map fields every time step but the issue was not fully resolved.

thanks alot,

James

[Lieven]

Hi Jim,

Did you manage to find a solution for this issue?

Cheers,

L

[carowjp]

Lieven,

Sorry for the delayed reply I was on a short vacation.

Unfortunately I have not been able to resolve this issue...but I haven't given up.

best regards,

Jim

[pixarzhang]

i am doing something with meshToMesh.interpolate now and encountered a problem. I used "meshToMeshInterp.interpolate(fieldTarget,fieldSour ce)"


,but when i wmake ,the terminal says

"meshandinit.H:71:65: error: no matching function for call to ‘Foam::meshToMesh::interpolate(Foam::volScalarFiel d&, Foam::volScalarField&)’
meshandinit.H:71:65: note: candidates are:
...."

it seems like there must be 4 candidates in my case.something like
{
GeometricField<Type, fvPatchField, volMesh>&,
const GeometricField<Type, fvPatchField, volMesh>&,
order,
const CombineOp& cop = eqOp<Type>()
}
when i looked into meshToMesh.H,it also has the CombineOp& cop

I try to add "meshToMesh::INTERPOLATE "in it, but I do not know "const CombineOp& cop = eqOp<Type>()"can set as what?





i am not very familiar with the source files. I also looked into the source files of mapField but do not find out what the cop is. do any body know where the CombineOp& cop comes from.


by the way,it seems that most people only used 2 candidates in interpolate like
"meshToMeshInterp.interpolate(fieldTarget,fieldSou rce)", but how can't i succeed.below is how i define the variables

fvMesh meshTarget
(
IOobject
(
fvMesh::defaultRegion,
runTimeTarget.timeName(),//runTimeTarget.constant()
runTimeTarget
)
);
meshToMesh meshToMeshInterp(mesh, meshTarget, patchMap_, cuttingPatches_);

volScalarField fieldSource
(
p,
mesh
);

volScalarField fieldTarget
(
p,
meshTarget
);

meshToMeshInterp.interpolate(fieldTarget,fieldSour ce);


any help would be appreciated.

[pixarzhang]

Originally Posted by pixarzhang

i am doing something with meshToMesh.interpolate now and encountered a problem. I used "meshToMeshInterp.interpolate(fieldTarget,fieldSour ce)"


,but when i wmake ,the terminal says

"meshandinit.H:71:65: error: no matching function for call to ‘Foam::meshToMesh::interpolate(Foam::volScalarFiel d&, Foam::volScalarField&)’
meshandinit.H:71:65: note: candidates are:
...."

it seems like there must be 4 candidates in my case.something like
{
GeometricField<Type, fvPatchField, volMesh>&,
const GeometricField<Type, fvPatchField, volMesh>&,
order,
const CombineOp& cop = eqOp<Type>()
}
when i looked into meshToMesh.H,it also has the CombineOp& cop

I try to add "meshToMesh::INTERPOLATE "in it, but I do not know "const CombineOp& cop = eqOp<Type>()"can set as what?





i am not very familiar with the source files. I also looked into the source files of mapField but do not find out what the cop is. do any body know where the CombineOp& cop comes from.


by the way,it seems that most people only used 2 candidates in interpolate like
"meshToMeshInterp.interpolate(fieldTarget,fieldSou rce)", but how can't i succeed.below is how i define the variables

fvMesh meshTarget
(
IOobject
(
fvMesh::defaultRegion,
runTimeTarget.timeName(),//runTimeTarget.constant()
runTimeTarget
)
);
meshToMesh meshToMeshInterp(mesh, meshTarget, patchMap_, cuttingPatches_);

volScalarField fieldSource
(
p,
mesh
);

volScalarField fieldTarget
(
p,
meshTarget
);

meshToMeshInterp.interpolate(fieldTarget,fieldSour ce);


any help would be appreciated.

I tried it out finally!!
the fourh candidate can set as " Foam::eqOp<scalar>() ",and the .C file can be compiled successfully.

despite i don't know where the class comes from ,i searched the src fold,but can't find any file named *eqOp* or related.

wish someone could tell me where to know what Foam::eqOp<scalar>() is .And all the questions before are also not resolved.

[wyldckat]

Quick answer: "Foam::eqOp<Type>()" is a special class with a single operator-method implemented, meant for helping with the process of combining data from several parallel processes.
By assigning a specific variable type in "Type", you specify which kind of data is being compared.

[pixarzhang]

Thank you for your explaination

but I am still not very clear.
I found a statement in meshToMeshinterplate.c which is "
const CombineOp& cop
...
cop(toF[celli],fromf[adr[celli]]) ",
I guess it give the value of fromf[adr[celli]] to toF[celli]?

do you know where the source file of the class is located, or where I can find it,I want to see more details if could.

[wyldckat]

Hi pixarzhang,

Which OpenFOAM version are you using? Because I'm looking at OpenFOAM 2.3.x and the code has changed considerably from your description. Therefore knowing which version of OpenFOAM you're using would make it a lot easier to explain this.

Best regards,
Bruno

[pixarzhang]

HI,wyldckat,
my OpenFOAM version is 2.2.0,I haven't updated it.
thanks again.

[wyldckat]

Hi pixarzhang,

I've finally managed to look into this. So what happens is this:

1. In OpenFOAM 2.1 it uses this line of code:
   

   toF[celli] = fromVf[adr[celli]];

2. In OpenFOAM 2.2, it uses this line of code:
   

   cop(toF[celli], fromVf[adr[celli]]);

3. The first one would simply state that the value "toF[celli]" on the current cell mesh "celli", should be the value that was on "fromVf[adr[celli]]" for the old cell address "adr[celli]".
4. The limitation of the first line of code is that it requires for both cell values (old and new locations) have to be in the same processor (aka sub-domain).
5. The second implementation, done in OpenFOAM 2.2, allows for the value to "fromVf[adr[celli]]" to come from another processor (aka sub-domain), which is why it is using "cop" as a combine operator, instead of standard equal operator.

This implementation apparently worked, but on OpenFOAM 2.3 they changed everything completely for this class "meshToMesh", apparently using a more robust and parallel way of doing its objectives.

Best regards,
Bruno

[pixarzhang]

thanks，that helps a lot

[Sylv]

To complete the answer of wyldckat, the meshtomesh class in 2.3.X is the only one (to my knowledge), which can permform the interpolation from a decomposted source mesh (srcMesh) to a decomposted target mesh (tgtMesh).

In my own solver, I interpolate the U field from a coarse source mesh (srcMesh) to a fine target mesh (tgtMesh). The tgtMesh represent only a small portion of the srcMesh. Here is the code, which does the job:

// a simple function
wordList addProcessorPatches
(
    const fvMesh& meshTarget,
    const wordList& cuttingPatches
)
{
    // Add the processor patches to the cutting list
    HashSet<word> cuttingPatchTable;
    forAll(cuttingPatches, i)
    {
        cuttingPatchTable.insert(cuttingPatches[i]);
    }

    const polyBoundaryMesh& pbm = meshTarget.boundaryMesh();

    forAll(pbm, patchI)
    {
        if (isA<processorPolyPatch>(pbm[patchI]))
        {
            const word& patchName = pbm[patchI].name();
            cuttingPatchTable.insert(patchName);
        }
    }

    return cuttingPatchTable.toc();
}

// generate the source mesh "srcMesh", the target mesh "tgtMesh",
// the source velocity field "U"  and the target velocity field 
// "Utgt" (fill it with vector::zero or some dummy value...)

A very long piece of code here. You should know how to do it....

// get patchMap (I read the hash table from a dict...)
HashTable<word> patchMap;
myDict.lookup("patchMap") >> patchMap;

// get the cutting patch. This the target patches which cut the source mesh
// I read this list from a dictionnary
wordList cuttingPatches;
myDict.lookup("cuttingPatches") >>  cuttingPatches;

// As I want to be able to run in parallel, the extra processor patches of the target
// mesh are also cutting patch, therefore, I need to update "cuttingPatch"
wordList allCuttingPatches;
allCuttingPatches = addProcessorPatches(tgtMesh, cuttingPatches)

//set a map method
meshToMesh::interpolationMethod mapMethod = meshToMesh::imCellVolumeWeight;

// create the meshtomesh interpolator object
meshToMesh interpolatorObj
(
            srcMesh,
            tgtMesh,
            mapMethod,
            patchMap, 
            allCuttingPatches
);

//do the interpolation between U of the source mesh and Utgt of the target mesh
interpolatorObj.mapSrcToTgt
        (
            U,
            eqOp<vector>(),
            Utgt
        );

I hope it will help

Marcel

[pixarzhang]

it does help.thanks~

[chriss85]

I'm now facing a similar task.
My solver uses a complete mesh that is then separated into different regions.
In the solver, I need to map a field between the full mesh and the separate regions.

I don't care about the boundary conditions of the fields as I can set them manually in the solver to zeroGradient. Do I need a patchMap and cuttingPatches then, or can I leave these lists empty (apart from the processor patches as you've shown before)?
If I leave them empty, do I need to manually assign a zeroGradient BC, or is it kept from the original field? If it is, a call to .correctBoundaryConditions() might be required.

I'm using OF 2.3.x, is the code you posted above valid for this version and parallel runs?

Also, for performance reasons, can we cache the mesh-to-mesh addressing? In my first tries the mapping has been quite slow.

[openfoam_aero]

Originally Posted by Sylv

To complete the answer of wyldckat, the meshtomesh class in 2.3.X is the only one (to my knowledge), which can permform the interpolation from a decomposted source mesh (srcMesh) to a decomposted target mesh (tgtMesh).

In my own solver, I interpolate the U field from a coarse source mesh (srcMesh) to a fine target mesh (tgtMesh). The tgtMesh represent only a small portion of the srcMesh. Here is the code, which does the job:

// a simple function
wordList addProcessorPatches
(
    const fvMesh& meshTarget,
    const wordList& cuttingPatches
)
{
    // Add the processor patches to the cutting list
    HashSet<word> cuttingPatchTable;
    forAll(cuttingPatches, i)
    {
        cuttingPatchTable.insert(cuttingPatches[i]);
    }

    const polyBoundaryMesh& pbm = meshTarget.boundaryMesh();

    forAll(pbm, patchI)
    {
        if (isA<processorPolyPatch>(pbm[patchI]))
        {
            const word& patchName = pbm[patchI].name();
            cuttingPatchTable.insert(patchName);
        }
    }

    return cuttingPatchTable.toc();
}

// generate the source mesh "srcMesh", the target mesh "tgtMesh",
// the source velocity field "U"  and the target velocity field 
// "Utgt" (fill it with vector::zero or some dummy value...)

A very long piece of code here. You should know how to do it....

// get patchMap (I read the hash table from a dict...)
HashTable<word> patchMap;
myDict.lookup("patchMap") >> patchMap;

// get the cutting patch. This the target patches which cut the source mesh
// I read this list from a dictionnary
wordList cuttingPatches;
myDict.lookup("cuttingPatches") >>  cuttingPatches;

// As I want to be able to run in parallel, the extra processor patches of the target
// mesh are also cutting patch, therefore, I need to update "cuttingPatch"
wordList allCuttingPatches;
allCuttingPatches = addProcessorPatches(tgtMesh, cuttingPatches)

//set a map method
meshToMesh::interpolationMethod mapMethod = meshToMesh::imCellVolumeWeight;

// create the meshtomesh interpolator object
meshToMesh interpolatorObj
(
            srcMesh,
            tgtMesh,
            mapMethod,
            patchMap, 
            allCuttingPatches
);

//do the interpolation between U of the source mesh and Utgt of the target mesh
interpolatorObj.mapSrcToTgt
        (
            U,
            eqOp<vector>(),
            Utgt
        );

I hope it will help

Marcel

Since it has been a few years and since this thread MASSIVELY helped me, I wish to contribute.

There are a few syntax changes that have taken place over the years. Particularly, the handling of the decomposition of the source and target meshes differently on the processors is improved thanks to the LOD and AABB methods.


I use Openfoam v1812 so please bear this in mind!

In light of that, we can change the code to

    meshToMesh InterpolatorObj
    (
    mesh,
    target,
    Foam::meshToMesh::interpolationMethod::imCellVolumeWeight,
    Foam::meshToMesh::procMapMethod::pmAABB,
    false
    );

    //do the interpolation between U of the source mesh and Utgt of the target mesh
    InterpolatorObj.mapSrcToTgt
    (
    U,
    plusEqOp<vector>(),
    U_target
    );

Here the only thing that has changed is that we add the procMapMethod pmAABB and false for interpolating all patches (I presume this is bnecause my case is not consistent so it would make sense to not interpolate all patches).

Hope this is helpful for someone who is working on this now!

I am calling this piece of code at the end of the simple loop (I am running a steady state case and I do not need the interpolations at every iteration). But as someone pointed out here, it would be nice to know how to have the cell addressing computed once and carry forward.