View Issue Details

IDProjectCategoryView StatusLast Update
0001347OpenFOAMBugpublic2015-03-05 17:53
Reporteruser965Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
PlatformLinuxOSScientific LinuxOS Version6.5 (Carbon)
Summary0001347: Computation of tetrahedra overlap volume fails on small grid cells
DescriptionThe method "tetOverlapVolume::tetTetOverlapVol(const tetPoints& tetA,
const tetPoints& tetB) const" always returns zero when the tetrahedra have a small volume. This makes the utility mapField fails on meshes with grid spacing < 1e-5. This method is called for instance by the interpolation method cellVolumeWeightMethod.
Steps To ReproduceUnpack the attached cases. The first case contains a 2x2 mesh with a field that has value 1. The second case contains a 3x3 mesh with exactly one overlapping grid cell. The latter contains a field with value "uniform 0". Run

$ ./Allrun
$ cat case2/0/testField

The utility mapFields with the method cellVolumeWeight gives the wrong result "uniform 0" (it actually does nothing), while the correct result is

internalField nonuniform List<scalar> 9(1 0 0 0 0 0 0 0 0);
Additional InformationThis bug is caused by the line 67 in tetOverlapVolume.C:

if ((tetA.tet().mag() < SMALL) || (tetB.tet().mag() < SMALL))
{
        return 0.0;
}

With double precision, SMALL = 1e-15. Therefore, when the cell size is 1e-5 or smaller, the resulting tetrahedrons always have a volume less than SMALL. This if-condition causes the overlap volume to be 0, which is very undesirable. Therefore, I suggest replacing SMALL by pow(SMALL, 3) (one could replace it by VSMALL, but that solution makes less sense to me).
TagsNo tags attached.

Activities

user965

2014-07-11 13:58

 

cases.tar (2,957 bytes)

user4

2015-03-05 09:33

 

Test-tetTetOverlap.C (4,273 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Application
    Test-tetTetOverlap

Description
    Overlap volume of two tets

\*---------------------------------------------------------------------------*/

#include "tetrahedron.H"
#include "OFstream.H"
#include "meshTools.H"

using namespace Foam;

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

void writeOBJ
(
    Ostream& os,
    label& vertI,
    const tetPoints& tet
)
{
    forAll(tet, fp)
    {
        meshTools::writeOBJ(os, tet[fp]);
    }
    os  << "l " << vertI+1 << ' ' << vertI+2 << nl
        << "l " << vertI+1 << ' ' << vertI+3 << nl
        << "l " << vertI+1 << ' ' << vertI+4 << nl
        << "l " << vertI+2 << ' ' << vertI+3 << nl
        << "l " << vertI+2 << ' ' << vertI+4 << nl
        << "l " << vertI+3 << ' ' << vertI+4 << nl;
    vertI += 4;
}


int main(int argc, char *argv[])
{
    Sout.precision(16);
    Pout.precision(16);

    tetPoints A
    (
        point(0, 0, 0),
        point(1, 0, 0),
        point(1, 1, 0),
        point(1, 1, 1)
    );
    const tetPointRef tetA = A.tet();

    tetPoints B
    (
        point(1, 0.1, 0.1),
        point(1.1, 0.1, 0.1),
        point(1.1, 1.1, 0.1),
        point(1.1, 1.1, 1.1)
    );
    const tetPointRef tetB = B.tet();


    tetPointRef::tetIntersectionList insideTets;
    label nInside = 0;
    tetPointRef::tetIntersectionList outsideTets;
    label nOutside = 0;

    tetA.tetOverlap
    (
        tetB,
        insideTets,
        nInside,
        outsideTets,
        nOutside
    );


    // Dump to file
    // ~~~~~~~~~~~~

    {
        OFstream str("tetA.obj");
        Info<< "Writing A to " << str.name() << endl;
        label vertI = 0;
        writeOBJ(str, vertI, A);
    }
    {
        OFstream str("tetB.obj");
        Info<< "Writing B to " << str.name() << endl;
        label vertI = 0;
        writeOBJ(str, vertI, B);
    }
    {
        OFstream str("inside.obj");
        Info<< "Writing parts of A inside B to " << str.name() << endl;
        label vertI = 0;
        for (label i = 0; i < nInside; ++i)
        {
            writeOBJ(str, vertI, insideTets[i]);
        }
    }
    {
        OFstream str("outside.obj");
        Info<< "Writing parts of A outside B to " << str.name() << endl;
        label vertI = 0;
        for (label i = 0; i < nOutside; ++i)
        {
            writeOBJ(str, vertI, outsideTets[i]);
        }
    }


    // Check
    // ~~~~~

    Info<< "Vol A:" << tetA.mag() << endl;

    scalar volInside = 0;
    for (label i = 0; i < nInside; ++i)
    {
        volInside += insideTets[i].tet().mag();
    }
    Info<< "Vol A inside B:" << volInside << endl;

    scalar volOutside = 0;
    for (label i = 0; i < nOutside; ++i)
    {
        volOutside += outsideTets[i].tet().mag();
    }
    Info<< "Vol A outside B:" << volOutside << endl;

    Info<< "Sum inside and outside:" << volInside+volOutside << endl;

    if (mag(volInside+volOutside-tetA.mag()) > SMALL)
    {
        FatalErrorIn("Test-tetetOverlap")
            << "Tet volumes do not sum up to input tet."
            << exit(FatalError);
    }

    return 0;
}


// ************************************************************************* //
Test-tetTetOverlap.C (4,273 bytes)   

user4

2015-03-05 09:38

  ~0003959

The test is quite strict. For two touching, regular, tets I see overlap volume:

    DP: 1.019927013371182e-47
    SP: 7.484330182379794e-21

However with making one of the tets have an aspect ratio of 1e-9 I get

    SP: 1.666666665789407e-10

Could you try

    sqr(SMALL)

instead of SMALL? This seems like a better compromise.

henry

2015-03-05 17:53

manager   ~0003969

commit fe29caafb0f0546c553391dba9377f08950284a6

Issue History

Date Modified Username Field Change
2014-07-11 13:58 user965 New Issue
2014-07-11 13:58 user965 File Added: cases.tar
2015-03-05 09:33 user4 File Added: Test-tetTetOverlap.C
2015-03-05 09:38 user4 Note Added: 0003959
2015-03-05 17:53 henry Note Added: 0003969
2015-03-05 17:53 henry Status new => resolved
2015-03-05 17:53 henry Resolution open => fixed
2015-03-05 17:53 henry Assigned To => henry