GGI and Topological changes - checkMesh fails

strakakl · June 18, 2012, 10:35

Hi foamers,

I am working to combine the GGI interface and topological changes of the mesh in OpenFOAM 1.6 ext. In the video you can see my test case, a cylinder where the right (or front) part is growing (translation) in the negative z-direction whereas the left (or back) part is rotating.

https://dl.dropbox.com/u/86131187/GgiTopoChange.avi

Calculation was done with the transientSimpleDyMFoam solver on a single processor, at first sight the results look ok. Post processing works fine with with the paraFoam utility.
But running checkMesh gives the following error

Code:

--> FOAM FATAL ERROR: 
Problem with patch-to zone addressing: some patch faces not found in interpolation zone

    From function void ggiPolyPatch::calcZoneAddressing() const
    in file meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C at line 103.

FOAM aborting

Aborted

After some debugging i found that there is something wrong with the addressing of patch faces of GGI patches (called "screw" and "chamber" in my test case). At time t=0.005 the face of the faceZone screw with the globalFaceID 129936 can not be found in the map of the localFaceIDs.

Code:

/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |
|  \\    /   O peration     | Version:  1.6-ext                               |
|   \\  /    A nd           | Web:      www.extend-project.de                 |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 1.6-ext-f3027b3161e4
Exec   : checkMesh
Date   : Jun 18 2012
Time   : 13:29:33
Host   : Brain
PID    : 8352
Case   : /home/ak113859/OpenFOAM/ak113859-1.6-ext/run/GgiTest/GgiWithTopoChanger
nProcs : 1
SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

--> FOAM Warning : 
    From function dlLibraryTable::open(const fileName& functionLibName)
    in file db/dlLibraryTable/dlLibraryTable.C at line 86
    could not load /home/ak113859/OpenFOAM/ak113859-1.6-ext/lib/linux64GccDPOpt/libipimdynamicFvMesh.so: undefined symbol: _ZTIN4Foam17topoChangerFvMeshE
Create polyMesh for time = 0

void polyMesh::initMesh() : initialising primitiveMesh
face zone name: screw
face zone index: 2
size: 2784
start: 121632
startLabel: 0
Initializing the GGI interpolator between master/shadow patches: outletScrew-7/inletChamber-10

    From function void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
    in file lnInclude/GGIInterpolationWeights.C at line 94
    Evaluation of GGI weighting factors:

    From function GGIInterpolation::findNonOverlappingFaces
    in file lnInclude/GGIInterpolationWeights.C at line 633
       : Found 0 non-overlapping faces for this GGI patch

    From function GGIInterpolation::findNonOverlappingFaces
    in file lnInclude/GGIInterpolationWeights.C at line 633
       : Found 0 non-overlapping faces for this GGI patch
  Largest slave weighting factor correction : 0.0105776 average: 0.000376867
  Largest master weighting factor correction: 0.0105228 average: 0.00035835
face zone name: chamber
face zone index: 0
size: 2784
start: 124416
startLabel: 0
polyMesh::globalData() const : Constructing parallelData from processor topology
polyMesh::readUpdateState polyMesh::readUpdate() : Updating mesh based on saved data.
Faces instance: old = "0" new = "0"
Points instance: old = "0" new = "0"
No change
Time = 0

void polyMesh::clearGeom() : clearing geometric data
void polyMesh::clearAddressing() : clearing topology
polyMesh::globalData() const : Constructing parallelData from processor topology
Mesh stats
    all points:           48490
    live points:           48490
    all faces:            130560
    live faces:            130560
    internal faces:   115968
    cells:            41088
    boundary patches: 7
    point zones:      0
    face zones:       3
    cell zones:       2

Overall number of cells of each type:
    hexahedra:     41088
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
   *Number of regions: 2
    The mesh has multiple regions which are not connected by any face.
  <<Writing region information to "0/cellToRegion"

Checking patch topology for multiply connected surfaces ...
    Patch               Faces    Points   Surface topology                  
    wallScrew-4         2784     2833     ok (non-closed singly connected)  
    wallBarrel-5        1536     1632     ok (non-closed singly connected)  
    inletScrew-6        1344     1440     ok (non-closed singly connected)  
    outletScrew-7       2784     2833     ok (non-closed singly connected)  
    inletChamber-10     2784     2833     ok (non-closed singly connected)  
    wallChamber-11      576      672      ok (non-closed singly connected)  
    front-12            2784     2833     ok (non-closed singly connected)  

Checking geometry...
    This is a 3-D mesh
    Overall domain bounding box (-0.00147418 -0.00147418 -0.00141191) (0.00147418 0.00147418 0.000669233)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Mesh (non-empty, non-wedge) dimensions 3
    Boundary openness (-2.33671e-18 -1.2622e-19 1.17466e-16) Threshold = 1e-06 OK.
    Max cell openness = 3.07599e-16 OK.
    Max aspect ratio = 6.80908 OK.
    Minumum face area = 4.34507e-10. Maximum face area = 9.66699e-09.  Face area magnitudes OK.
    Min volume = 2.96127e-14. Max volume = 4.2747e-13.  Total volume = 9.5907e-09.  Cell volumes OK.
    Mesh non-orthogonality Max: 34.5777 average: 3.28214 Threshold = 70
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.80492 OK.

Mesh OK.

polyMesh::readUpdateState polyMesh::readUpdate() : Updating mesh based on saved data.
Faces instance: old = "0" new = "0.005"
Points instance: old = "0" new = "0.005"
Topological change
void polyMesh::clearGeom() : clearing geometric data
void polyMesh::clearAddressing() : clearing topology
void polyMesh::setInstance(const fileName& inst) : Resetting file instance to "0.005"
void polyMesh::setMotionWriteOpt(IOobject::writeOption) Setting motion writeOpt to 0
void polyMesh::setTopoWriteOpt(IOobject::writeOption) Setting topo writeOpt to 0
void polyMesh::initMesh() : initialising primitiveMesh
face zone name: screw
face zone index: 2
size: 2784
start: 129936
startLabel: -1


--> FOAM FATAL ERROR: 
Problem with patch-to zone addressing: some patch faces not found in interpolation zone

    From function void ggiPolyPatch::calcZoneAddressing() const
    in file meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C at line 103.

FOAM aborting

Aborted

faceZone at t=0.005 reads

Code:

...
screw
{
    type faceZone;
faceLabels      List<label> 
2784
(
129936
129937
...
132719
)
;
flipMap         List<bool> 2784{1};
}
)

Further investigations show that none of the faces of "screw" can be found in the map. Taking a look on the map of local face indices it seems that this map contains at time step t=0.005 the informations of time step t=0, because the local indices reaches from 121632 to 124415.

faceZone at time step 0

Code:

...
screw
{
    type faceZone;
faceLabels      List<label> 
2784
(
121632
121633
...
124415
)
;
flipMap         List<bool> 2784{1};
}
)

It looks like that i have to update this map at each new time step but at the moment i don't have any idea to do that.

Someone out there who can shed some light on this issue?

Best regards,

Klaus

strakakl · June 19, 2012, 12:14

some update to my post above

running checkMesh utility for a single time step works perfect

Code:

checkMesh -time 0

or

Code:

checkMesh -time 0.005

results both in Mesh OK. Indeed it seems that there is some update of the faceZones missing.

Today i tried to find out what's going wrong. In the file polyMeshIO.C the line

Code:

boundary_.calcGeometry();

is the last one executed before the fatal error. Taking a closer look at this file one can see that point-, face-, and cell zones are updated afterwards.

Code:

...
        // Calculate topology for the patches (processor-processor comms etc.)
        boundary_.updateMesh();
        
        // Calculate the geometry for the patches (transformation tensors etc.)
        boundary_.calcGeometry();
        
        // Derived info
        bounds_ = boundBox(allPoints_);
        geometricD_ = Vector<label>::zero;
        solutionD_ = Vector<label>::zero;

        // Zones
        pointZoneMesh newPointZones
        (
            IOobject
            (
                "pointZones",
                facesInst,
                meshSubDir,
                *this,
                IOobject::READ_IF_PRESENT,
                IOobject::NO_WRITE,
                false
            ),
            *this
        );

        label oldSize = pointZones_.size();

        if (newPointZones.size() <= pointZones_.size())
        {
            pointZones_.setSize(newPointZones.size());
        }

        // Reset existing ones
        forAll (pointZones_, czI)
        {
            pointZones_[czI] = newPointZones[czI];
        }

        // Extend with extra ones
        pointZones_.setSize(newPointZones.size());

        for (label czI = oldSize; czI < newPointZones.size(); czI++)
        {
            pointZones_.set(czI, newPointZones[czI].clone(pointZones_));
        }

        pointZones_.setSize(newPointZones.size());
        forAll (pointZones_, pzI)
        {
            pointZones_[pzI] = newPointZones[pzI];
        }

        faceZoneMesh newFaceZones
        (
            IOobject
            (
                "faceZones",
                facesInst,
                meshSubDir,
                *this,
                IOobject::READ_IF_PRESENT,
                IOobject::NO_WRITE,
                false
            ),
            *this
        );

        oldSize = faceZones_.size();

        if (newFaceZones.size() <= faceZones_.size())
        {
            faceZones_.setSize(newFaceZones.size());
        }

        // Reset existing ones
        forAll (faceZones_, fzI)
        {
            faceZones_[fzI].resetAddressing
            (
                newFaceZones[fzI],
                newFaceZones[fzI].flipMap()
            );
        }

        // Extend with extra ones
        faceZones_.setSize(newFaceZones.size());

        for (label fzI = oldSize; fzI < newFaceZones.size(); fzI++)
        {
            faceZones_.set(fzI, newFaceZones[fzI].clone(faceZones_));
        }


        cellZoneMesh newCellZones
        (
            IOobject
            (
                "cellZones",
                facesInst,
                meshSubDir,
                *this,
                IOobject::READ_IF_PRESENT,
                IOobject::NO_WRITE,
                false
            ),
            *this
        );

        oldSize = cellZones_.size();

        if (newCellZones.size() <= cellZones_.size())
        {
            cellZones_.setSize(newCellZones.size());
        }

        // Reset existing ones
        forAll (cellZones_, czI)
        {
            cellZones_[czI] = newCellZones[czI];
        }

        // Extend with extra ones
        cellZones_.setSize(newCellZones.size());

        for (label czI = oldSize; czI < newCellZones.size(); czI++)
        {
            cellZones_.set(czI, newCellZones[czI].clone(cellZones_));
        }
...

Moving the boundary_.* stuff behind the update of the zones

Code:

...
        // Zones
        pointZoneMesh newPointZones
        (
            IOobject
            (
                "pointZones",
                facesInst,
                meshSubDir,
                *this,
                IOobject::READ_IF_PRESENT,
                IOobject::NO_WRITE,
                false
            ),
            *this
        );

        label oldSize = pointZones_.size();

        if (newPointZones.size() <= pointZones_.size())
        {
            pointZones_.setSize(newPointZones.size());
        }

        // Reset existing ones
        forAll (pointZones_, czI)
        {
            pointZones_[czI] = newPointZones[czI];
        }

        // Extend with extra ones
        pointZones_.setSize(newPointZones.size());

        for (label czI = oldSize; czI < newPointZones.size(); czI++)
        {
            pointZones_.set(czI, newPointZones[czI].clone(pointZones_));
        }

        pointZones_.setSize(newPointZones.size());
        forAll (pointZones_, pzI)
        {
            pointZones_[pzI] = newPointZones[pzI];
        }

        faceZoneMesh newFaceZones
        (
            IOobject
            (
                "faceZones",
                facesInst,
                meshSubDir,
                *this,
                IOobject::READ_IF_PRESENT,
                IOobject::NO_WRITE,
                false
            ),
            *this
        );

        oldSize = faceZones_.size();

        if (newFaceZones.size() <= faceZones_.size())
        {
            faceZones_.setSize(newFaceZones.size());
        }

        // Reset existing ones
        forAll (faceZones_, fzI)
        {
            faceZones_[fzI].resetAddressing
            (
                newFaceZones[fzI],
                newFaceZones[fzI].flipMap()
            );
        }

        // Extend with extra ones
        faceZones_.setSize(newFaceZones.size());

        for (label fzI = oldSize; fzI < newFaceZones.size(); fzI++)
        {
            faceZones_.set(fzI, newFaceZones[fzI].clone(faceZones_));
        }


        cellZoneMesh newCellZones
        (
            IOobject
            (
                "cellZones",
                facesInst,
                meshSubDir,
                *this,
                IOobject::READ_IF_PRESENT,
                IOobject::NO_WRITE,
                false
            ),
            *this
        );

        oldSize = cellZones_.size();

        if (newCellZones.size() <= cellZones_.size())
        {
            cellZones_.setSize(newCellZones.size());
        }

        // Reset existing ones
        forAll (cellZones_, czI)
        {
            cellZones_[czI] = newCellZones[czI];
        }

        // Extend with extra ones
        cellZones_.setSize(newCellZones.size());

        for (label czI = oldSize; czI < newCellZones.size(); czI++)
        {
            cellZones_.set(czI, newCellZones[czI].clone(cellZones_));
        }

         // Calculate topology for the patches (processor-processor comms etc.)
         boundary_.updateMesh();
         
         // Calculate the geometry for the patches (transformation tensors etc.)
         boundary_.calcGeometry();
         
         // Derived info
         bounds_ = boundBox(allPoints_);
         geometricD_ = Vector<label>::zero;
         solutionD_ = Vector<label>::zero;
...

solves the problem, checkMesh now results in Mesh OK

Question to the experts: is it allowed to do that or are there any side effects which i should be aware of? Comments would be fine.

Best regards,

Klaus

deepsterblue · June 19, 2012, 12:37

Hmm... looks like the zones were written to disk with addressing prior to topo-changes. This may not be an actual problem, since you say that the solution looks fine, but an issue where inconsistent addressing is written to disk. Still needs to be fixed though...

strakakl · June 20, 2012, 06:30

Quote:

Hmm... looks like the zones were written to disk with addressing prior to topo-changes.

i think the zones are written to disk correctly, the faceZone files are modified after a topo-change.

t=0;

Code:

... screw 
{     
    type faceZone; 
    faceLabels      List<label>  2784 
    (
         121632 
         121633 
         ...
         124415 
     ) 
     ;
     flipMap         List<bool> 2784{1};
 } 
....

-> topo-Change t=0.005

Code:

... screw
    {
    type faceZone; 
    faceLabels      List<label>  
    2784 
    ( 
    129936
    129937 
    ...
    132719 
    )
    ;
    flipMap         List<bool> 2784{1};
    } 
...

As you can see the face IDs are different after the topoChange.
The problem arises when the mesh is read from disk.

Best,

Klaus

Hgholami · November 10, 2019, 06:48

Dear foamer
Did you involve with this error?

Quote:

FatalErrorIn Reconstructed cell centres already calculated occurs.

it is become from
Quote:
// Create neighbouring face centres using interpolation
if (master())
{
const label shadowID = shadowIndex();

// Get the transformed and interpolated shadow face cell centers
reconFaceCellCentresPtr_ =
new vectorField
(
interpolate
(
boundaryMesh()[shadowID].faceCellCentres()
- boundaryMesh()[shadowID].faceCentres()
)
+ faceCentres()
);
}
else
{
FatalErrorIn("void ggiPolyPatch::calcReconFaceCellCentres() const")
<< "Attempting to create reconFaceCellCentres on a shadow"
<< abort(FatalError);
}
Do you know, why this error appears?
when I use cyclicGgi for top and bottom of rectangular domain, the problem appear when running simulation.

June 19, 2012, 12:37		#3
deepsterblue Senior Member Sandeep Menon Join Date: Mar 2009 Location: Amherst, MA Posts: 403 Rep Power: 25	Hmm... looks like the zones were written to disk with addressing prior to topo-changes. This may not be an actual problem, since you say that the solution looks fine, but an issue where inconsistent addressing is written to disk. Still needs to be fixed though... __________________ Sandeep Menon University of Massachusetts Amherst https://github.com/smenon

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Create a GeometricField of a given type on given patch?	philippose	OpenFOAM Programming & Development	4	August 12, 2013 13:41
engine Simulations with GGI instead of Sliding Interface	abminternet	OpenFOAM	3	April 30, 2013 04:41
GGI together with topological changes...	kalle	OpenFOAM Running, Solving & CFD	6	May 7, 2012 14:50
GGI modeling	shavkat	OpenFOAM Running, Solving & CFD	2	December 22, 2010 05:39