CFD Online Discussion Forums - Mapping of pointFields with topology changes

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Mapping of pointFields with topology changes (https://www.cfd-online.com/Forums/openfoam-programming-development/106469-mapping-pointfields-topology-changes.html)

Mapping of pointFields with topology changes

Hello all,

This question is related to the thread http://www.cfd-online.com/Forums/ope...ent-2-1-x.html, but since this is more general it might not be found later on by interested parties.

When a mesh undergoes a topology change, the updateFields() function is called for the mesh classes in a way that should map all fields in the object registry. In testing combinations of mesh topology/movement operations, I've found that only volume and surface fields are mapped by default. The following code calls the functions that should do the mapping:

Code:

pointMesh motionUpdate(*this);

motionUpdate.updateMesh(map);

, where *this is a mesh class pointer. However, executing this gives a message such as the following for all pointFields

Code:

Not mapping point<Type>Field <fieldName> since originating mesh differs from that of mapper.

Does anyone know how to properly do this mapping, or if it is even possible in the current 2.1.x release?

Thanks!

You shouldn't have to worry about this. Point fields should be mapped automatically in the polyMesh::updateMesh(const mapPolyMesh& mpm) member function. Take a look at polyMeshUpdate.C. The pointMesh should automatically map all related pointFields (see pointMesh.C and the pointMeshMapper for details).

Thanks for the reply Sandeep; I assume you are referring to this code snippet:

Code:

// pointMesh

if (thisDb().foundObject<pointMesh>(pointMesh::typeName))

{

     const_cast<pointMesh&>

     (

         thisDb().lookupObject<pointMesh>

         (

             pointMesh::typeName

         )

     ).updateMesh(mpm);

}

which definitely looks like it should do the mapping for pointfields. That is why the message that the pointField isn't being updated so puzzling. Besides not registering the pointField with the objectRegistry (which doesn't make sense as a volScalarField created inside the motionSolver IS mapped automatically), can you think of a reason why the pointMesh wouldn't get mapped? The message says that the mapping is refused as the originating mesh differs from the mapper; I tried placing the pointMesh update after the updateMesh for refinement and it still fails.

I've added the files required to run the case I've been testing it against (I hadn't included a solver; my bad), as well as commented the lines in dynamicMotionSolverRefineFvMesh::refine and unrefine that try to do the mapping. Thanks for the help!

It appears that MapGeometricFields does a pointer comparison of the mesh returned by the mapper and the mesh defined for the field. Perhaps you could check through gdb to see which one of the two is inconsistent?

Print out the addresses of both the polyMesh and pointMesh objects, and it should become clear. If there's an inconsistency somewhere, then this is surely a bug - possibly due to a polyMesh reference being returned instead of a pointMesh, or something similar.

Ah.. I get what might be wrong.. You're instantiating a pointMesh within the updateMesh member function within dynamicRefineFvMesh, which destroys itself when the function goes out of scope.

Obviously, any pointField instantiated prior to this call is using a pointMesh registered elsewhere, which is why it doesn't get mapped. You _have_ to use the polyMesh update functionality - there's simple no other way.

Hope this helps.

That makes sense; I was attempting the manual method to provide what is seemingly missing from the code. What is odd is that right after the refinement there is a call to updateMesh(map), which following the chain up should end up in the pointField mapping and execute, even if it only returns that "unable to map" message. I think this may be a bug: will post a report with the relevant files and see what happens.

EDIT: Here is the bug report, if anyone finds anything relevant be sure to share it with the developers: http://www.openfoam.org/mantisbt/view.php?id=638

Dear Marco and Sandeep

I am trying to add Mesh refinment to my code using multiDirRefinement.
The code is doing the refinment on newMesh which is at the begining the same as mesh .
then with changemesh the mesh is changed according to newMesh.

Code:

multiDirRefinement multiRef(newMesh, refCells, refineDict);

        polyTopoChange meshMod(newMesh);

       autoPtr<mapPolyMesh> morphMapPtr = meshMod.changeMesh(mesh, false, true);

       const mapPolyMesh& morphMap = morphMapPtr();

But I dont know how to map the feilds which have been calculated on the mesh ( not the newMesh) to the mesh which is changed.

Best
Mahdi

Hello Mahdi,

After calling changemesh, you are returned a mapPolyMesh that needs to be applied to the fields. There should be a member function for your mesh called "updateMesh" that takes as its argument the mapPolyMesh. Have a look at the post earlier up where there is a code snippet for updating the pointmesh. Most meshes that hold data have this member function.

Dear Marco

Thanks for your answer. I have tried your suggestion.

Code:

multiDirRefinement multiRef(newMesh, refCells, refineDict);

 newMesh.write();

 newMesh.readUpdate();

 polyTopoChange meshMod(newMesh);

 autoPtr<mapPolyMesh> morphMapPtr = meshMod.changeMesh(mesh, false);

 mesh.write();

 const mapPolyMesh& morphMap = morphMapPtr();

//mesh.readUpdate();

mesh.updateMesh(morphMap);

But it gives me a following error :

Code:

#0  Foam::error::printStack(Foam::Ostream&) in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"

#1  Foam::sigFpe::sigFpeHandler(int) in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"

#2   in "/lib/libc.so.6"

#3  Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"

#4  

 in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/myicoFoam"

#5  

 in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/myicoFoam"

#6  __libc_start_main in "/lib/libc.so.6"

#7  

 in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/myicoFoam"

Floating point exception

Best
Mahdi

Does this happen at the solver step? Check to make sure the fields were actually mapped and then we can see what the problem may be.

Dear Marco

Thanks for the reply.

Yes its exactly happening at the solver step . field are not mapped as you guessed. but the mesh is refined.

Best
Mahdi

Quote:

Originally Posted by mturcios777 (Post 390650)

Does this happen at the solver step? Check to make sure the fields were actually mapped and then we can see what the problem may be.

Dear Marco

I have changed the position of the refinement and mapping to after the solution of equations.
The data are now mapped. but still there is a problem. I have checked the values for p. some of them are almost NAN and some of them zero. the mapping of the data has problem. and because of that I am receiving the following error.

Code:

#0  Foam::error::printStack(Foam::Ostream&) in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"

#1  Foam::sigFpe::sigFpeHandler(int) in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"

#2   in "/lib/libc.so.6"

#3  Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"

#4  

 in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/myicoFoam"

#5  

 in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/myicoFoam"

#6  __libc_start_main in "/lib/libc.so.6"

#7  

 in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/myicoFoam"

Floating point exception

Best
Mahdi

If p is nan/0 in locations then there is something wrong with the mapPolyMesh that was created. I seem to recall that when using refinement that some of the fields are reconstructed (like flux and such). Have a look at dynamicRefineFvMesh.C (particularly the refine/unrefine functions) to see what else needs to be done to properly map all the fields.

On that subject, is it only p that is mapped improperly, or are other fields still problematic?

Good luck,
Marco

UPDATE: Its been a while, but the bug reported above has been fixed. Coupled with some changes to the pEqn in the engine solvers, layerAdditionRemoval can now be used! You will need to make a new engineMesh class that has the polyTopoChangers included. This is truly awesome! Still trying to get it to work with dynamicFvMeshes, but there appears to be an issue with the thermophysical properties coupling. More updates on other threads as we go.

Dear Marco

I have checked that the problem is somehow both at p and U. I have doubt as you mention that my map maybe wrong. the mesh modifer that i am using is multiDirRefinment which is using the undoablemeshcutter when doing the refinment on specified direction. the thing is that when doing the refinment with dynamicRefineFvMesh it is using the hexRef8 meshcutter. the diffrence is that the hexRef8 meshcutter will use the poloytopochange (meshmod) as input for mesh refinment. but in the multiDirRefinment it uses the mesh.
because of that i use to equal mesh at the beginging to generate the map. doing the refinment on one and the changing the the other according to poloytopochange of the other and creating the map.

Code:

    multiDirRefinement multiRef(newMesh, refCells, refineDict);

     newMesh.write();

    multiDirRefinement multiRef(newMesh, refCells, refineDict);

     newMesh.write();

     polyTopoChange meshMod(newMesh);

     autoPtr<mapPolyMesh> morphMapPtr = meshMod.changeMesh(mesh,false);

     const mapPolyMesh& morphMap = morphMapPtr();

     mesh.updateMesh(morphMap);  

    mesh.write();

the second point is that I also uses the meshtomesh interpolation between my two different meshes. doing the refinement just on newMesh, interpolate the data to the newMesh from Mesh for storing the values. then doing the refinement on mesh and interpolating the data to mesh using mesh to mesh interpolation. untill here p and U are mapped correctly . but the problem is interpolating the surfacescaler as meshtomesh interpolation cannot handle that. so i try to do construct the phi by interpolating velocities. but it gives a error . infact no data is created for phi !!!

/////////////////
MyIcoFoam

Code:



#include "fvCFD.H"

#include "multiDirRefinement.H" 

#include "undoableMeshCutter.H" 

#include "polyTopoChange.H"

#include "errorEstimate.H"

#include "resError.H"

#include "evaluateError.H"

#include "errorDrivenRefinement.H"

#include "meshToMesh.H"

#include "GeometricField.H"

#include "IOobjectList.H"

#include "mapPolyMesh.H"

#include "polyMesh.H"

#include "dynamicRefineFvMesh.H"

#include "addToRunTimeSelectionTable.H"

#include "surfaceFields.H"

#include "pointFields.H"



static const scalar errTol = 1E-5;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



int main(int argc, char *argv[])

{

    #include "setRootCase.H"

    #include "createTime.H"

    #include "createmainMesh.H"

   #include "createNewMesh.H"

 //   #include "createDynamicFvMesh.H"

    #include "createFields.H"

   #include "createFieldsn.H"

    #include "initContinuityErrs.H"

    Info<< "\nStarting time loop\n" << endl;

    

while (runTime.run())

    {



        scalar  realTime=0.0;

        Info<< "Time = " << runTime.timeName() << nl << endl;



        #include "readPISOControls.H"

        #include "CourantNo.H"



if(runTime.write()) 

{

//////////////////////////////////////////////////////////////////////////

///  Error Calculation 



   #include "IcoFoamErrorCalc.H"



/////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////

/// Mark Cells For refinment



   #include "Markcellforref.H"



/////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////

///  Refinment



 #include "RefinmentwithoutMap.H"

/////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////

///  Mapping 

  #include "MaptomeshNew.H"

  #include "MapfromMeshnew.H"

} 



runTime++;



    fvVectorMatrix UEqn

        (

            fvm::ddt(U)

          + fvm::div(phi, U)

          - fvm::laplacian(nu, U)

        );



        solve(UEqn == -fvc::grad(p));



        // --- PISO loop



        for (int corr=0; corr<nCorr; corr++)

        {

            volScalarField rUA = 1.0/UEqn.A();



            U = rUA*UEqn.H();

            phi = (fvc::interpolate(U) & mesh.Sf());

               + fvc::ddtPhiCorr(rUA, U, phi);



            adjustPhi(phi, U, p);



            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)

            {

                fvScalarMatrix pEqn

                (

                    fvm::laplacian(rUA, p) == fvc::div(phi)

                );



                pEqn.setReference(pRefCell, pRefValue);

                pEqn.solve();



                if (nonOrth == nNonOrthCorr)

                {

                    phi -= pEqn.flux();

                }

            }

            #include "continuityErrs.H"



            U -= rUA*fvc::grad(p);

            U.correctBoundaryConditions();

        }

runTime.write();



        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

            << "  ClockTime = " << runTime.elapsedClockTime() << " s"

            << nl << endl;

    }



    Info<< "End\n" << endl;



    return 0;

}

////////////////
interpolating from mesh to meshNew

Code:

         



meshToMesh  meshToMeshInterp(mesh, newMesh);

        const fvMesh& meshSource = meshToMeshInterp.fromMesh();

        const fvMesh& meshTarget = meshToMeshInterp.toMesh();

       

        volScalarField pnn

        (

        IOobject

        (

            "pnn",

            runTime.timeName(),

            meshTarget,

            IOobject::READ_IF_PRESENT,

            IOobject::AUTO_WRITE

        ),

        meshTarget,

        dimensionedScalar("zero", dimensionSet(0,2,-2,0,0,0,0), scalar(0.0))

        );





        volVectorField Unn

        (

        IOobject

        (

            "Unn",

            runTime.timeName(),

            meshTarget,

            IOobject::READ_IF_PRESENT,

            IOobject::AUTO_WRITE

        ),

        meshTarget,

        dimensionedVector("zero", dimensionSet(0,1,-1,0,0,0,0), vector::zero)

        );



  // Interpolate field



Info << "***************interpolate from mesh to meshNEw ***********************" << endl;

Info << "***********  interpolate p to pnn ***************" << endl;

       meshToMeshInterp.interpolate(pnn,p,meshToMesh::MAP);



         pnn.write();

Info << "***********  interpolate U to Unn ***************" << endl;

       meshToMeshInterp.interpolate(Unn,U,meshToMesh::MAP);



         Unn.write();

Info << "***********  interpolation Finish***************" << endl;

Info << "***********  mesh change according to meshmod of newMesh***************" << endl;

     polyTopoChange meshMod(newMesh);

     autoPtr<mapPolyMesh> morphMapPtr = meshMod.changeMesh(mesh, false);

     mesh.write();

//////////////////////////
interpolating from meshNew to mesh

Code:

   meshToMesh  meshToMeshInterpb(newMesh, mesh);

    const fvMesh& meshSourceb = meshToMeshInterpb.fromMesh();

    const fvMesh& meshTargetb = meshToMeshInterpb.toMesh();



      GeometricField<scalar, fvPatchField, volMesh> pfieldSourceb

        (

            pnn,

            meshSourceb

        );



    p=dimensionedScalar("zero", dimensionSet(0,2,-2,0,0,0,0), scalar(0.0));

    p.write();



 IOobject fieldTargetIOobjectb

        (

            "p",

            meshTargetb.time().timeName(),

            mesh,

            IOobject::READ_IF_PRESENT,

            IOobject::AUTO_WRITE

        );

            // Read fieldTarget

            GeometricField<scalar, fvPatchField, volMesh> fieldTargetb

            (

                fieldTargetIOobjectb,

                meshTargetb

            ); 



      GeometricField<vector, fvPatchField, volMesh> ufieldSourceb

        (

            Unn,

            meshSourceb

        );



    U=dimensionedVector("zero", dimensionSet(0,1,-1,0,0,0,0), vector::zero);

    U.write();



       IOobject fieldTargetIOobjectub

        (

            "U",

            meshTargetb.time().timeName(),

            meshTargetb,

            IOobject::MUST_READ,

            IOobject::AUTO_WRITE

        );



            // Read fieldTarget

            GeometricField<vector, fvPatchField, volMesh> fieldTargetub

            (

                fieldTargetIOobjectub,

                meshTargetb

            ); 

 Info << "***************interpolate from meshNEw to mesh ***********************" << endl;



// Interpolate field

Info << "***********  interpolate pnn to p ***************" << endl;

           meshToMeshInterpb.interpolate(fieldTargetb,pnn,meshToMesh::MAP);

            fieldTargetb.write();

Info << "***********  interpolate Unn to U ***************" << endl;

           meshToMeshInterpb.interpolate(fieldTargetub,Unn,meshToMesh::MAP);

            fieldTargetub.write();



Info << "***********  interpolate phi on surface ***************" << endl;



       phi=dimensionedScalar("zero", dimensionSet(0,3,-1,0,0,0,0), scalar(0.0));

       phi.write();

       phi=(dimensionedVector("zero", dimensionSet(0,1,-1,0,0,0,0), vector::zero) & mesh.Sf());

       phi.write();

       phi=(fvc::interpolate(U) & mesh.Sf());

       phi.write();



Info << "***********" << endl;

///////////////////////////////
Error

Code:

#0  Foam::error::printStack(Foam::Ostream&) in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"

#1  Foam::sigSegv::sigSegvHandler(int) in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"

#2   in "/lib/libc.so.6"

#3  Foam::surfaceInterpolationScheme<Foam::Vector<double> >::interpolate(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::tmp<Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> > const&) in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/myicoFoam"

#4  Foam::surfaceInterpolationScheme<Foam::Vector<double> >::interpolate(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&) const in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libfiniteVolume.so"

#5  

 at myicoFoam.C:0

#6  

 in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/myicoFoam"

#7  __libc_start_main in "/lib/libc.so.6"

#8  

 in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/myicoFoam"

Segmentation fault

if i move this part of code (refinment and mapping ) to the end of code after runTime.write() , some value are mapped for phi . the values are wrong and even there is some value for boundaries. the error that this time is occurring is

Code:



Foam::error::printStack(Foam::Ostream&) in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"

#1  Foam::sigFpe::sigFpeHandler(int) in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"

#2   in "/lib/libc.so.6"

#3  Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"

#4  void Foam::divide<Foam::fvsPatchField, Foam::surfaceMesh>(Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh>&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/myicoFoam"

#5  

 in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/myicoFoam"

#6  

 in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/myicoFoam"

#7  __libc_start_main in "/lib/libc.so.6"

#8  

 in "/home/mahdi/OpenFOAM/OpenFOAM-1.7.1/applications/bin/linux64GccDPOpt/myicoFoam"

Floating point exception

I will be thankful if i have your help.

Best
Mahdi

Dear Marco

I have also added #include "CorrectFluxinterpolated.H" to my code according to dynamicRefineFvMesh.C . this code is supposed to check and correct the phi while it is mapped.!!!!?
But now just it is correcting the boundary condition and make the value of phi zero. moreover there is still a NAN in one cell surface in phi.

Best
Mahdi

Dear Mahdi

I've got the same error info as you said when foam does flux update if my guess is correct after UpdateMesh(map).

May I know whether you have fixed the problem or not and how? thanks in advance!

Best Regards

sxh

Dear Foamers,

I am struggling with a similar problem for the past few days.

My aim is to implement a dynamic mesh into the dsmc Solver, which is able to adapt the mesh to the current flow. Here it is necessary to refine the mesh in specified areas. I have chosen a similar procedure like it is described above. My code looks like:

Code:



    fvMesh newMesh      

    (

      Foam::IOobject

      (

          Foam::fvMesh::defaultRegion,

          runTime.timeName(),

          runTime,

          Foam::IOobject::MUST_READ,

          Foam::IOobject::AUTO_WRITE

      )

    );

    

    multiDirRefinement multiRef(newMesh,refineCells(),refineDict);

    

    polyTopoChange changer(newMesh);

    

    autoPtr<mapPolyMesh> map = changer.changeMesh(mesh_, false);

    const mapPolyMesh& m = map();

    

    autoMap(m);

I am starting with creating a new mesh which then is refined. (using multiDirRefinement and polyTopoChange like described above)
Then I use the mapPolyMesh to do the mapping. The first part of the function autoMap() looks like:

Code:

    typedef typename  ParcelType::trackingData tdType;

    tdType td(*this);

    Cloud<ParcelType>::template autoMap<tdType>(td, mapper);

    

    mesh_.updateMesh(mapper);

The error occurs when the mapping of the internalField's is done inside updateMesh(). Following error message:

Code:

--> FOAM FATAL ERROR: 

index 2816 out of range 0 ... 2815



    From function void Foam::UList<T>::checkIndex(Foam::label) const [with T = double; Foam::label = int]

    in file /home/ip56/OpenFOAM/OpenFOAM-dev/src/OpenFOAM/lnInclude/UListI.H at line 109.



FOAM aborting

The old mesh has a cell number of 2816 and the new Mesh 5917. I am now wondering why the mapper has a cell storage of size 2816.

Am I missing some steps? I am happy for any suggestions.

Best regards