diff --git a/ThirdParty/rpmBuild/SOURCES/ParMGridGen-1.0.patch b/ThirdParty/rpmBuild/SOURCES/ParMGridGen-1.0.patch index 1bf93c7ad..394e0b346 100644 --- a/ThirdParty/rpmBuild/SOURCES/ParMGridGen-1.0.patch +++ b/ThirdParty/rpmBuild/SOURCES/ParMGridGen-1.0.patch @@ -1,6 +1,31 @@ +diff -ruN ParMGridGen-1.0_orig/Makefile ParMGridGen-1.0/Makefile +--- ParMGridGen-1.0_orig/Makefile 2017-04-04 15:02:44.020713666 +0200 ++++ ParMGridGen-1.0/Makefile 2017-04-04 15:21:48.582647336 +0200 +@@ -1,16 +1,21 @@ + default: ++ (mkdir bin) + (cd MGridGen ; make) + + serial: ++ (mkdir bin) + (cd MGridGen ; make) + + parallel: ++ (mkdir bin) + (cd MGridGen ; make) + (cd ParMGridGen ; make) + clean: ++ (mkdir bin) + (cd MGridGen ; make clean) + (cd ParMGridGen ; make clean ) + + realclean: ++ (mkdir bin) + (cd MGridGen ; make realclean ) + (cd ParMGridGen ; make realclean ) diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in ---- ParMGridGen-1.0_orig/Makefile.in 2001-12-04 16:30:33.000000000 -0800 -+++ ParMGridGen-1.0/Makefile.in 2013-08-22 20:07:33.491171127 -0700 +--- ParMGridGen-1.0_orig/Makefile.in 2017-04-04 15:02:44.012713543 +0200 ++++ ParMGridGen-1.0/Makefile.in 2017-04-04 15:06:00.159742074 +0200 @@ -1,6 +1,6 @@ #-------------------------------------------------------------------------- # Which make to use @@ -18,6 +43,15 @@ diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in # Which loader to use LD = cc +@@ -22,7 +22,7 @@ + LDOPTIONS = -O3 + + # Where to put the executable +-BINDIR = ../.. ++BINDIR = ../../bin + + # Additional libraries + DMALLOCDIR = /usr/local @@ -33,22 +33,25 @@ # In which directories to look for any additional libraries diff --git a/ThirdParty/rpmBuild/SOURCES/ParMGridGen-1.0.patch_darwin b/ThirdParty/rpmBuild/SOURCES/ParMGridGen-1.0.patch_darwin index 29a9db97c..2f793fb40 100644 --- a/ThirdParty/rpmBuild/SOURCES/ParMGridGen-1.0.patch_darwin +++ b/ThirdParty/rpmBuild/SOURCES/ParMGridGen-1.0.patch_darwin @@ -79,8 +79,8 @@ diff -ruN ParMGridGen-1.0_orig/MGridGen/Programs/Makefile ParMGridGen-1.0/MGridG ifeq ($(ddmalloc),yes) DEBUGFLAGS := $(DEBUGFLAGS) -DDMALLOC -DDEBUG diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in ---- ParMGridGen-1.0_orig/Makefile.in 2011-12-24 13:54:44.000000000 -0500 -+++ ParMGridGen-1.0/Makefile.in 2011-12-24 13:49:26.000000000 -0500 +--- ParMGridGen-1.0_orig/Makefile.in 2001-12-05 01:30:33.000000000 +0100 ++++ ParMGridGen-1.0/Makefile.in 2017-04-04 15:36:04.695980033 +0200 @@ -1,6 +1,6 @@ #-------------------------------------------------------------------------- # Which make to use @@ -98,6 +98,15 @@ diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in # Which loader to use LD = cc +@@ -22,7 +22,7 @@ + LDOPTIONS = -O3 + + # Where to put the executable +-BINDIR = ../.. ++BINDIR = ../../bin + + # Additional libraries + DMALLOCDIR = /usr/local @@ -33,22 +33,25 @@ # In which directories to look for any additional libraries @@ -129,3 +138,28 @@ diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in #-------------------------------------------------------------------------- # +diff -ruN ParMGridGen-1.0_orig/Makefile ParMGridGen-1.0/Makefile +--- ParMGridGen-1.0_orig/Makefile 2001-11-09 00:41:22.000000000 +0100 ++++ ParMGridGen-1.0/Makefile 2017-04-04 14:51:04.033914737 +0200 +@@ -1,16 +1,21 @@ + default: ++ (mkdir bin) + (cd MGridGen ; make) + + serial: ++ (mkdir bin) + (cd MGridGen ; make) + + parallel: ++ (mkdir bin) + (cd MGridGen ; make) + (cd ParMGridGen ; make) + clean: ++ (mkdir bin) + (cd MGridGen ; make clean) + (cd ParMGridGen ; make clean ) + + realclean: ++ (mkdir bin) + (cd MGridGen ; make realclean ) + (cd ParMGridGen ; make realclean ) diff --git a/ThirdParty/rpmBuild/SOURCES/ParaView-4.4.0.patch_darwin b/ThirdParty/rpmBuild/SOURCES/ParaView-4.4.0.patch_darwin new file mode 100644 index 000000000..3d1d8ce4b --- /dev/null +++ b/ThirdParty/rpmBuild/SOURCES/ParaView-4.4.0.patch_darwin @@ -0,0 +1,12 @@ +diff -ruN ParaView-4.3.1_orig/Applications/ParaView-4.3.1_extra_install_Darwin.cmake ParaView-4.3.1/Applications/ParaView-4.3.1_extra_install_Darwin.cmake +--- ParaView-4.3.1_orig/Applications/ParaView-4.3.1_extra_install_Darwin.cmake 1969-12-31 19:00:00.000000000 -0500 ++++ ParaView-4.3.1/Applications/ParaView-4.3.1_extra_install_Darwin.cmake 2013-10-02 19:00:00.000000000 -0400 +@@ -0,0 +1,8 @@ ++# ++# Additional install rules for Mac OS X platforms ++# ++INSTALL (DIRECTORY buildObj/bin/paraview.app ++ DESTINATION ${CMAKE_INSTALL_PREFIX}/bin ++ USE_SOURCE_PERMISSIONS ++ COMPONENT Runtime) ++ diff --git a/ThirdParty/rpmBuild/SOURCES/libccmio-2.6.1.patch_0 b/ThirdParty/rpmBuild/SOURCES/libccmio-2.6.1.patch_0 index dbbcb47aa..4a82e7fcf 100644 --- a/ThirdParty/rpmBuild/SOURCES/libccmio-2.6.1.patch_0 +++ b/ThirdParty/rpmBuild/SOURCES/libccmio-2.6.1.patch_0 @@ -1,13 +1,13 @@ diff -ruN libccmio-2.6.1_orig/config/config.gnu.to.star libccmio-2.6.1/config/config.gnu.to.star ---- libccmio-2.6.1_orig/config/config.gnu.to.star 2007-12-17 18:32:03.000000000 -0500 -+++ libccmio-2.6.1/config/config.gnu.to.star 2010-12-18 14:50:35.000000000 -0500 +--- libccmio-2.6.1_orig/config/config.gnu.to.star 2007-12-18 00:32:03.000000000 +0100 ++++ libccmio-2.6.1/config/config.gnu.to.star 2017-04-04 13:41:35.000000000 +0200 @@ -1,4 +1,4 @@ -#!/bin/sh +#!/bin/bash # $Id: config.gnu.to.star,v 1.4 2006/06/05 21:12:16 geoffp Exp $ -@@ -34,6 +34,12 @@ +@@ -34,6 +34,13 @@ x86_64-unknown-linux-gnu-null) echo linux64_2.4-x86-glibc_2.2.5 ;; ppc64-unknown-linux-gnu-null) echo linux64_2.6-pwr4-glibc_2.3.3 ;; i386-apple-darwin8-null) echo i386-apple-darwin8 ;; @@ -17,6 +17,7 @@ diff -ruN libccmio-2.6.1_orig/config/config.gnu.to.star libccmio-2.6.1/config/co + i386-apple-darwin13-null) echo i386-apple-darwin13 ;; + i386-apple-darwin14-null) echo i386-apple-darwin14 ;; + i386-apple-darwin15-null) echo i386-apple-darwin15 ;; ++ i386-apple-darwin16-null) echo i386-apple-darwin16 ;; *) echo unknown ;; esac @@ -30,15 +31,15 @@ diff -ruN libccmio-2.6.1_orig/config/config.guess libccmio-2.6.1/config/config.g # Copyright (C) 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001 # Free Software Foundation, Inc. diff -ruN libccmio-2.6.1_orig/config/config.system libccmio-2.6.1/config/config.system ---- libccmio-2.6.1_orig/config/config.system 2008-02-25 22:07:16.000000000 -0500 -+++ libccmio-2.6.1/config/config.system 2010-12-18 14:51:34.000000000 -0500 +--- libccmio-2.6.1_orig/config/config.system 2008-02-26 04:07:16.000000000 +0100 ++++ libccmio-2.6.1/config/config.system 2017-04-04 13:42:15.000000000 +0200 @@ -1,4 +1,4 @@ -#! /bin/sh +#! /bin/bash # $Id: config.system,v 1.2 2005/09/29 22:19:19 geoffp Exp $ -@@ -87,6 +87,27 @@ +@@ -87,6 +87,30 @@ i386-apple-darwin8.11.1) echo i386-apple-darwin8 ;; @@ -62,6 +63,9 @@ diff -ruN libccmio-2.6.1_orig/config/config.system libccmio-2.6.1/config/config. + + i386-apple-darwin15.* ) + echo i386-apple-darwin15 ;; ++ ++ i386-apple-darwin16.* ) ++ echo i386-apple-darwin16 ;; + *) echo unknown diff --git a/ThirdParty/rpmBuild/SPECS/libccmio-2.6.1.spec b/ThirdParty/rpmBuild/SPECS/libccmio-2.6.1.spec index a8c47a8da..b62d21d0f 100644 --- a/ThirdParty/rpmBuild/SPECS/libccmio-2.6.1.spec +++ b/ThirdParty/rpmBuild/SPECS/libccmio-2.6.1.spec @@ -105,6 +105,7 @@ Patch0: libccmio-2.6.1.patch_0 [ ! -d config/i386-apple-darwin13 ] && cp -r config/i386-apple-darwin8 config/i386-apple-darwin13 [ ! -d config/i386-apple-darwin14 ] && cp -r config/i386-apple-darwin8 config/i386-apple-darwin14 [ ! -d config/i386-apple-darwin15 ] && cp -r config/i386-apple-darwin8 config/i386-apple-darwin15 + [ ! -d config/i386-apple-darwin16 ] && cp -r config/i386-apple-darwin8 config/i386-apple-darwin16 %endif # Warning: # 1: The name of the ADF library will be renamed to libadf_ccmio since this diff --git a/applications/utilities/mesh/conversion/star4ToFoam/cellZoneToCellTableId/cellZoneToCellTableId.C b/applications/utilities/mesh/conversion/star4ToFoam/cellZoneToCellTableId/cellZoneToCellTableId.C index 8de3f0332..ed7bf3baa 100644 --- a/applications/utilities/mesh/conversion/star4ToFoam/cellZoneToCellTableId/cellZoneToCellTableId.C +++ b/applications/utilities/mesh/conversion/star4ToFoam/cellZoneToCellTableId/cellZoneToCellTableId.C @@ -27,6 +27,7 @@ Description \*---------------------------------------------------------------------------*/ #include "argList.H" +#include "objectRegistry.H" #include "foamTime.H" #include "volFields.H" #include "OFstream.H" diff --git a/applications/utilities/mesh/conversion/star4ToFoam/foamMeshToStar/foamMeshToStar.C b/applications/utilities/mesh/conversion/star4ToFoam/foamMeshToStar/foamMeshToStar.C index 565204c3a..fb25291ef 100644 --- a/applications/utilities/mesh/conversion/star4ToFoam/foamMeshToStar/foamMeshToStar.C +++ b/applications/utilities/mesh/conversion/star4ToFoam/foamMeshToStar/foamMeshToStar.C @@ -35,6 +35,7 @@ Description \*---------------------------------------------------------------------------*/ #include "argList.H" +#include "objectRegistry.H" #include "foamTime.H" #include "volFields.H" #include "cellModeller.H" diff --git a/etc/bashrc b/etc/bashrc index c0cf57639..47ae6e22d 100755 --- a/etc/bashrc +++ b/etc/bashrc @@ -160,10 +160,11 @@ done : ${WM_OSTYPE:=POSIX}; export WM_OSTYPE -# Compiler: set to Gcc or Icc (for Intel's icc) +# Compiler: set to Gcc, Icc (for Intel's icc) or Clang # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ : ${WM_COMPILER:=Gcc}; export WM_COMPILER #: ${WM_COMPILER:=Icc}; export WM_COMPILER +#: ${WM_COMPILER:=Clang}; export WM_COMPILER export WM_COMPILER_ARCH= export WM_COMPILER_LIB_ARCH= @@ -379,7 +380,7 @@ Darwin) echo "Using Macports binaries" fi - export WM_USE_MACPORT=1 + export WM_USE_MACPORT=0 export WM_BASE_COMPILER=`echo $WM_COMPILER | tr -d "[:digit:]"` export WM_MACPORT_MPI_VERSION=`echo $WM_COMPILER | tr "[:upper:]" "[:lower:]"` export WM_MACPORT_VERSION=`echo $WM_MACPORT_MPI_VERSION | tr -d "[:alpha:]" | sed -e "s/\(.\)\(.\)/\1\.\2/"` @@ -448,8 +449,14 @@ Darwin) export WM_FC="gfortran-mp-$WM_MACPORT_VERSION" elif [ "$WM_BASE_COMPILER" == "Clang" ] then - export WM_CC="clang-mp-$WM_MACPORT_VERSION" - export WM_CXX="clang++-mp-$WM_MACPORT_VERSION" + if [ "$WM_USE_MACPORT" == "1" ] + then + export WM_CC="clang-mp-$WM_MACPORT_VERSION" + export WM_CXX="clang++-mp-$WM_MACPORT_VERSION" + else + export WM_CC="clang" + export WM_CXX="clang++" + fi # Seems like there is no Fortran-frontend for LLVM at thistime elif [ "$WM_BASE_COMPILER" == "Dragonegg" ] then diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C index fe4628f75..4d818a408 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C @@ -163,8 +163,7 @@ tmp > ggiFvPatchField::patchNeighbourField() const // Get shadow face-cells and assemble shadow field // This is a patchInternalField of neighbour but access is inconvenient. - // Assemble by hand. - // HJ, 27/Sep/2011 + // Assemble by hand. HJ, 27/Sep/2011 const unallocLabelList& sfc = ggiPatch_.shadow().faceCells(); Field sField(sfc.size()); @@ -180,13 +179,15 @@ tmp > ggiFvPatchField::patchNeighbourField() const if (ggiPatch_.bridgeOverlap()) { // Symmetry treatment used for overlap - vectorField nHat = this->patch().nf(); + const vectorField nHat = this->patch().nf(); // Use mirrored neighbour field for interpolation // HJ, 21/Jan/2009 - Field bridgeField = + const Field bridgeField = transform(I - 2.0*sqr(nHat), this->patchInternalField()); + // Note: bridging now takes into account fully uncovered and partially + // covered faces. VV, 18/Oct/2017. ggiPatch_.bridge(bridgeField, pnf); } @@ -211,18 +212,8 @@ void ggiFvPatchField::initEvaluate + (1.0 - this->patch().weights())*this->patchNeighbourField() ); - if (ggiPatch_.bridgeOverlap()) - { - // Symmetry treatment used for overlap - vectorField nHat = this->patch().nf(); - - Field pif = this->patchInternalField(); - - Field bridgeField = - 0.5*(pif + transform(I - 2.0*sqr(nHat), pif)); - - ggiPatch_.bridge(bridgeField, pf); - } + // Note: bridging and correction of partially overlapping faces taken into + // account in patchNeighbourField(). VV, 16/Oct/2017. Field::operator=(pf); } @@ -266,6 +257,21 @@ void ggiFvPatchField::initInterfaceMatrixUpdate scalarField pnf = ggiPatch_.interpolate(sField); + if (ggiPatch_.bridgeOverlap()) + { + // Note: will not work properly for types with rank > 0 (everything + // above scalar) if the symmetry plane is not aligned with one of the + // coordinate axes. VV, 18/Oct/2017. + const scalarField bridgeField = + transform + ( + I - 2.0*sqr(this->patch().nf()), + ggiPatch_.patchInternalField(psiInternal) + ); + + ggiPatch_.bridge(bridgeField, pnf); + } + // Multiply the field by coefficients and add into the result const unallocLabelList& fc = ggiPatch_.faceCells(); @@ -327,6 +333,21 @@ void ggiFvPatchField::initInterfaceMatrixUpdate Field pnf = ggiPatch_.interpolate(sField); + if (ggiPatch_.bridgeOverlap()) + { + // Note: will not work properly for types with rank > 0 (everything + // above scalar) if the symmetry plane is not aligned with one of the + // coordinate axes. VV, 18/Oct/2017. + const Field bridgeField = + transform + ( + I - 2.0*sqr(this->patch().nf()), + ggiPatch_.patchInternalField(psiInternal) + ); + + ggiPatch_.bridge(bridgeField, pnf); + } + // Multiply neighbour field with coeffs and re-use pnf for result // of multiplication multiply(pnf, coeffs, pnf); diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C index 3b1addbde..344b568c2 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C @@ -62,14 +62,14 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const // HJ, 2/Aug/2007 if (ggiPolyPatch_.master()) { - vectorField n = nf(); + const vectorField n = nf(); // Note: mag in the dot-product. // For all valid meshes, the non-orthogonality will be less than // 90 deg and the dot-product will be positive. For invalid // meshes (d & s <= 0), this will stabilise the calculation // but the result will be poor. HJ, 24/Aug/2011 - scalarField nfc = + const scalarField nfc = mag(n & (ggiPolyPatch_.reconFaceCellCentres() - Cf())); w = nfc/(mag(n & (Cf() - Cn())) + nfc); @@ -78,7 +78,9 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const { // Set overlap weights to 0.5 and use mirrored neighbour field // for interpolation. HJ, 21/Jan/2009 - bridge(scalarField(size(), 0.5), w); + const scalarField bridgedField(size(), 0.5); + + bridge(bridgedField, w); } } else @@ -87,7 +89,7 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const scalarField masterWeights(shadow().size()); shadow().makeWeights(masterWeights); - scalarField oneMinusW = 1 - masterWeights; + const scalarField oneMinusW = 1 - masterWeights; w = interpolate(oneMinusW); @@ -95,7 +97,9 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const { // Set overlap weights to 0.5 and use mirrored neighbour field // for interpolation. HJ, 21/Jan/2009 - bridge(scalarField(size(), 0.5), w); + const scalarField bridgedField(size(), 0.5); + + bridge(bridgedField, w); } } } @@ -107,16 +111,12 @@ void Foam::ggiFvPatch::makeDeltaCoeffs(scalarField& dc) const if (ggiPolyPatch_.master()) { // Stabilised form for bad meshes. HJ, 24/Aug/2011 - vectorField d = delta(); + const vectorField d = delta(); dc = 1.0/max(nf() & d, 0.05*mag(d)); - if (bridgeOverlap()) - { - scalarField bridgeDeltas = nf() & fvPatch::delta(); - - bridge(bridgeDeltas, dc); - } + // Note: no need to bridge the overlap since delta already takes it into + // account. VV, 18/Oct/2017. } else { @@ -126,7 +126,11 @@ void Foam::ggiFvPatch::makeDeltaCoeffs(scalarField& dc) const if (bridgeOverlap()) { - scalarField bridgeDeltas = nf() & fvPatch::delta(); + // Note: double the deltaCoeffs because this is symmetry treatment + // and fvPatch::deltaCoeffs() is cell to face inverse distance, + // while we need cell to "symmetry neighbour cell" distance. + // VV, 18/Oct/2017. + const scalarField bridgeDeltas = 2.0*fvPatch::deltaCoeffs(); bridge(bridgeDeltas, dc); } @@ -143,8 +147,8 @@ void Foam::ggiFvPatch::makeCorrVecs(vectorField& cv) const // Calculate correction vectors on coupled patches const scalarField& patchDeltaCoeffs = deltaCoeffs(); - vectorField patchDeltas = delta(); - vectorField n = nf(); + const vectorField patchDeltas = delta(); + const vectorField n = nf(); // If non-orthogonality is over 90 deg, kill correction vector // HJ, 6/Jan/2011 @@ -161,7 +165,10 @@ Foam::tmp Foam::ggiFvPatch::delta() const if (bridgeOverlap()) { - vectorField bridgeDeltas = Cf() - Cn(); + // Note: double the deltas because this is symmetry treatment and + // fvPatch::delta() is cell to face distance, while we need cell to + // "symmetry neighbour cell" distance. VV, 18/Oct/2017. + const vectorField bridgeDeltas = 2.0*fvPatch::delta(); bridge(bridgeDeltas, tDelta()); } @@ -177,7 +184,10 @@ Foam::tmp Foam::ggiFvPatch::delta() const if (bridgeOverlap()) { - vectorField bridgeDeltas = Cf() - Cn(); + // Note: double the deltas because this is symmetry treatment and + // fvPatch::delta() is cell to face distance, while we need cell to + // "symmetry neighbour cell" distance. VV, 18/Oct/2017. + const vectorField bridgeDeltas = 2.0*fvPatch::delta(); bridge(bridgeDeltas, tDelta()); } diff --git a/src/foam/interpolations/GGIInterpolation/GGIInterpolate.C b/src/foam/interpolations/GGIInterpolation/GGIInterpolate.C index b5fa94572..2e15382a6 100644 --- a/src/foam/interpolations/GGIInterpolation/GGIInterpolate.C +++ b/src/foam/interpolations/GGIInterpolation/GGIInterpolate.C @@ -97,13 +97,26 @@ void GGIInterpolation::bridge ( const Field& bridgeField, Field& ff, - const labelList& addr + const labelList& addr, + const labelList& partiallyCoveredAddr, + const scalarField& coveredFractions ) { + // Fully uncovered faces forAll (addr, faceI) { ff[addr[faceI]] = bridgeField[addr[faceI]]; } + + // Loop through partially covered faces and correct them. Note the + // operator+= since we assume that the interpolation part is carried out + // before bridging (see e.g. ggiFvPatchField::patchNeighbourField()) using + // weights that do not sum up to 1 + forAll (partiallyCoveredAddr, pcfI) + { + ff[partiallyCoveredAddr[pcfI]] += + coveredFractions[pcfI]*bridgeField[partiallyCoveredAddr[pcfI]]; + } } @@ -114,7 +127,9 @@ void GGIInterpolation::maskedBridge const Field& bridgeField, Field& ff, const labelList& mask, - const labelList& uncoveredFaces + const labelList& uncoveredFaces, + const labelList& partiallyCoveredAddr, + const scalarField& coveredFractions ) { // Note: tricky algorithm @@ -131,7 +146,7 @@ void GGIInterpolation::maskedBridge const label faceI = uncoveredFaces[uncoI]; // Search through the mask - for (; maskAddrI < mask.size(); maskAddrI++) + for (; maskAddrI < mask.size(); ++maskAddrI) { if (faceI == mask[maskAddrI]) { @@ -147,7 +162,38 @@ void GGIInterpolation::maskedBridge // Go one back and check for next uncovered face if (maskAddrI > 0) { - maskAddrI--; + --maskAddrI; + } + + break; + } + } + } + + // Reset maskAddrI + maskAddrI = 0; + + forAll (partiallyCoveredAddr, pcfI) + { + // Pick partially covered face + const label faceI = partiallyCoveredAddr[pcfI]; + + for (; maskAddrI < mask.size(); ++maskAddrI) + { + if (faceI == mask[maskAddrI]) + { + // Found masked partially covered face + ff[maskAddrI] += coveredFractions[pcfI]*bridgeField[maskAddrI]; + + break; + } + else if (mask[maskAddrI] > faceI) + { + // Gone beyond my index: my face is not present in the mask + // Go one back and check for next uncovered face + if (maskAddrI > 0) + { + --maskAddrI; } break; @@ -511,7 +557,8 @@ void GGIInterpolation::bridgeMaster if ( bridgeField.size() != masterPatch_.size() - || ff.size() != masterPatch_.size()) + || ff.size() != masterPatch_.size() + ) { FatalErrorIn ( @@ -527,7 +574,14 @@ void GGIInterpolation::bridgeMaster << abort(FatalError); } - bridge(bridgeField, ff, uncoveredMasterFaces()); + bridge + ( + bridgeField, + ff, + uncoveredMasterFaces(), + partiallyCoveredMasterFaces(), + masterFaceCoveredFractions() + ); } @@ -551,7 +605,7 @@ void GGIInterpolation::maskedBridgeMaster " Field& ff,\n" " const labelList& mask\n" ") const" - ) << "given field does not correspond to patch. Patch size: " + ) << "given field does not correspond to patch. Patch (mask) size: " << masterPatch_.size() << " bridge field size: " << bridgeField.size() << " field size: " << ff.size() @@ -559,7 +613,15 @@ void GGIInterpolation::maskedBridgeMaster << abort(FatalError); } - maskedBridge(bridgeField, ff, mask, uncoveredMasterFaces()); + maskedBridge + ( + bridgeField, + ff, + mask, + uncoveredMasterFaces(), + partiallyCoveredMasterFaces(), + masterFaceCoveredFractions() + ); } @@ -592,7 +654,14 @@ void GGIInterpolation::bridgeSlave << abort(FatalError); } - bridge(bridgeField, ff, uncoveredSlaveFaces()); + bridge + ( + bridgeField, + ff, + uncoveredSlaveFaces(), + partiallyCoveredSlaveFaces(), + slaveFaceCoveredFractions() + ); } @@ -616,7 +685,7 @@ void GGIInterpolation::maskedBridgeSlave " Field& ff\n," " const labelList& mask\n" ") const" - ) << "given field does not correspond to patch. Patch size: " + ) << "given field does not correspond to patch. Patch (mask) size: " << slavePatch_.size() << " bridge field size: " << bridgeField.size() << " field size: " << ff.size() @@ -624,12 +693,18 @@ void GGIInterpolation::maskedBridgeSlave << abort(FatalError); } - maskedBridge(bridgeField, ff, mask, uncoveredSlaveFaces()); + maskedBridge + ( + bridgeField, + ff, + mask, + uncoveredSlaveFaces(), + partiallyCoveredSlaveFaces(), + slaveFaceCoveredFractions() + ); } - - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/foam/interpolations/GGIInterpolation/GGIInterpolation.C b/src/foam/interpolations/GGIInterpolation/GGIInterpolation.C index 3bf5d55bd..8204b0bb7 100644 --- a/src/foam/interpolations/GGIInterpolation/GGIInterpolation.C +++ b/src/foam/interpolations/GGIInterpolation/GGIInterpolation.C @@ -104,7 +104,11 @@ void GGIInterpolation::clearOut() deleteDemandDrivenData(slaveWeightsPtr_); deleteDemandDrivenData(uncoveredMasterAddrPtr_); + deleteDemandDrivenData(partiallyCoveredMasterAddrPtr_); + deleteDemandDrivenData(masterFaceCoveredFractionsPtr_); deleteDemandDrivenData(uncoveredSlaveAddrPtr_); + deleteDemandDrivenData(partiallyCoveredSlaveAddrPtr_); + deleteDemandDrivenData(slaveFaceCoveredFractionsPtr_); deleteDemandDrivenData(masterPointAddressingPtr_); deleteDemandDrivenData(masterPointWeightsPtr_); @@ -155,7 +159,11 @@ GGIInterpolation::GGIInterpolation slavePointWeightsPtr_(NULL), slavePointDistancePtr_(NULL), uncoveredMasterAddrPtr_(NULL), - uncoveredSlaveAddrPtr_(NULL) + partiallyCoveredMasterAddrPtr_(NULL), + masterFaceCoveredFractionsPtr_(NULL), + uncoveredSlaveAddrPtr_(NULL), + partiallyCoveredSlaveAddrPtr_(NULL), + slaveFaceCoveredFractionsPtr_(NULL) { // Check size of transform. They should be equal to slave patch size // if the transform is not constant @@ -256,6 +264,32 @@ GGIInterpolation::uncoveredMasterFaces() const } +template +const labelList& +GGIInterpolation::partiallyCoveredMasterFaces() const +{ + if (!partiallyCoveredMasterAddrPtr_) + { + calcAddressing(); + } + + return *partiallyCoveredMasterAddrPtr_; +} + + +template +const scalarField& +GGIInterpolation::masterFaceCoveredFractions() const +{ + if (!masterFaceCoveredFractionsPtr_) + { + calcAddressing(); + } + + return *masterFaceCoveredFractionsPtr_; +} + + template const labelList& GGIInterpolation::uncoveredSlaveFaces() const @@ -269,6 +303,32 @@ GGIInterpolation::uncoveredSlaveFaces() const } +template +const labelList& +GGIInterpolation::partiallyCoveredSlaveFaces() const +{ + if (!partiallyCoveredSlaveAddrPtr_) + { + calcAddressing(); + } + + return *partiallyCoveredSlaveAddrPtr_; +} + + +template +const scalarField& +GGIInterpolation::slaveFaceCoveredFractions() const +{ + if (!slaveFaceCoveredFractionsPtr_) + { + calcAddressing(); + } + + return *slaveFaceCoveredFractionsPtr_; +} + + template bool GGIInterpolation::movePoints ( diff --git a/src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H b/src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H index 8a5e9e598..48dcf2f1c 100644 --- a/src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H +++ b/src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H @@ -219,9 +219,30 @@ class GGIInterpolation //- List of uncovered master faces mutable labelList* uncoveredMasterAddrPtr_; + //- List of partially covered master faces. Contains faces with + // weights between masterNonOverlapFaceTol_ and + // uncoveredFaceAreaTol_. Used in bridging + mutable labelList* partiallyCoveredMasterAddrPtr_; + + //- List of face fractions for partially covered master faces. Range + // anywhere from 1 to uncoveredFaceAreaTol_ (e.g. 0.2 means + // that the face area is covered by 20% on the other side) + mutable scalarField* masterFaceCoveredFractionsPtr_; + //- List of uncovered slave faces mutable labelList* uncoveredSlaveAddrPtr_; + //- List of partially covered slave faces. Contains faces with + // weights between masterNonOverlapFaceTol_ and + // uncoveredFaceAreaTol_. Used in bridging + mutable labelList* partiallyCoveredSlaveAddrPtr_; + + //- List of face fractions for partially covered slave faces. Range + // anywhere from 1 to uncoveredFaceAreaTol_ (e.g. 0.2 means + // that the face area is covered by 20% on the other side) + mutable scalarField* slaveFaceCoveredFractionsPtr_; + + // Private static data //- Facet area error tolerance @@ -233,6 +254,9 @@ class GGIInterpolation //- Facet bound box extension factor static const debug::tolerancesSwitch faceBoundBoxExtendSpanFraction_; + //- Tolerance for uncovered face areas + static const debug::tolerancesSwitch uncoveredFaceAreaTol_; + //- The next 3 attributes are parameters controlling the creation // of an octree search engine for the GGI facets neighbourhood // determination. @@ -417,6 +441,14 @@ class GGIInterpolation const scalar& nonOverlapFaceTol ) const; + //- Calculate partially covered faces + void calcPartiallyCoveredFaces + ( + const scalarListList& patchWeights, + const scalar& nonOverlapFaceTo, + const bool isMaster + ) const; + //- Clear all geometry and addressing void clearOut(); @@ -451,7 +483,9 @@ class GGIInterpolation ( const Field& bridgeField, Field& ff, - const labelList& addr + const labelList& addr, + const labelList& partiallyCoveredAddr, + const scalarField& coveredFractions ); //- Bridge uncovered faces given addressing @@ -462,7 +496,9 @@ class GGIInterpolation const Field& bridgeField, Field& ff, const labelList& mask, - const labelList& uncoveredFaces + const labelList& uncoveredFaces, + const labelList& partiallyCoveredAddr, + const scalarField& coveredFractions ); //- Is a transform required? @@ -533,9 +569,21 @@ public: //- Return reference to the master list of non-overlap faces const labelList& uncoveredMasterFaces() const; + //- Return reference to master list of partially covered faces + const labelList& partiallyCoveredMasterFaces() const; + + //- Return covered fractions of partially covered master faces + const scalarField& masterFaceCoveredFractions() const; + //- Return reference to the slave list of non-overlap faces const labelList& uncoveredSlaveFaces() const; + //- Return reference to slave list of partially covered faces + const labelList& partiallyCoveredSlaveFaces() const; + + //- Return covered fractions of partially covered slave faces + const scalarField& slaveFaceCoveredFractions() const; + //- Return reference to point addressing const List& masterPointAddr() const; @@ -554,6 +602,7 @@ public: //- Return distance to intersection for patch points const scalarField& slavePointDistanceToIntersection() const; + // Interpolation functions //- Interpolate from master to slave diff --git a/src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C b/src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C index 48845d702..4b3225498 100644 --- a/src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C +++ b/src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C @@ -68,6 +68,16 @@ GGIInterpolation::featureCosTol_ ); +template +const Foam::debug::tolerancesSwitch +GGIInterpolation::uncoveredFaceAreaTol_ +( + "GGIUncoveredFaceAreaTol", + 0.999, + "Fraction of face area mismatch (sum of weights) to consider a face " + "as uncovered, i.e. not to rescale weights." +); + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // template @@ -80,7 +90,11 @@ void GGIInterpolation::calcAddressing() const || slaveAddrPtr_ || slaveWeightsPtr_ || uncoveredMasterAddrPtr_ + || partiallyCoveredMasterAddrPtr_ + || masterFaceCoveredFractionsPtr_ || uncoveredSlaveAddrPtr_ + || partiallyCoveredSlaveAddrPtr_ + || slaveFaceCoveredFractionsPtr_ ) { FatalErrorIn @@ -431,7 +445,7 @@ void GGIInterpolation::calcAddressing() const ); // Fix: previously checked for VSMALL. - // HJ, 19/Sep2/106 + // HJ, 19/Sep/2016 if (intersectionArea > intersectionTestArea) { // We compute the GGI weights based on this @@ -555,7 +569,7 @@ void GGIInterpolation::calcAddressing() const const labelList& curMa = ma[mfI]; const scalarList& cursmaW = smaW[mfI]; - // Current master face indes + // Current master face index const label faceMaster = mfI; forAll (curMa, mAI) @@ -601,6 +615,29 @@ void GGIInterpolation::calcAddressing() const findNonOverlappingFaces(saW, slaveNonOverlapFaceTol_) ); + // Calculate master and slave partially covered addressing + + // Note: This function allocates: + // 1. partiallyCoveredMasterAddrPtr_ + // 2. masterFaceCoveredFractionsPtr_ + calcPartiallyCoveredFaces + ( + maW, + masterNonOverlapFaceTol_, + true // This is master + ); + + // Note: this function allocates: + // 1. partiallyCoveredSlaveAddrPtr_ + // 2. slaveFaceCoveredFractionsPtr_ + calcPartiallyCoveredFaces + ( + saW, + slaveNonOverlapFaceTol_, + false // This is not master + ); + + // Rescaling the weighting factors so they will sum up to 1.0 // See the comment for the method ::rescaleWeightingFactors() for // more information. By default, we always rescale. But for some @@ -608,7 +645,8 @@ void GGIInterpolation::calcAddressing() const // then we need the brute values, so no rescaling in that // case. Hence the little flag rescaleGGIWeightingFactors_ - // Not parallelised. HJ, 27/Apr/2016 + // Not parallelised. HJ, 27/Apr/2016. Rescaling is not performed for + // partially overlapping faces for their correct treatment. VV, 16/Oct/2017. if (rescaleGGIWeightingFactors_) { rescaleWeightingFactors(); @@ -656,6 +694,35 @@ void GGIInterpolation::rescaleWeightingFactors() const scalar sumMWC = 0; scalar curMWC = 0; + // Note: do not rescale weighting factors for partially covered faces + if (!partiallyCoveredMasterAddrPtr_ || !partiallyCoveredSlaveAddrPtr_) + { + FatalErrorIn + ( + "void GGIInterpolation::" + "rescaleWeightingFactors() const" + ) << "Master or slave partially covered faces are not calculated." + << abort(FatalError); + } + + const labelList& partiallyCoveredMasterFaces = + *partiallyCoveredMasterAddrPtr_; + const labelList& partiallyCoveredSlaveFaces = + *partiallyCoveredSlaveAddrPtr_; + + // Create a mask for partially covered master/slave faces + boolList masterPCMask(maW.size(), false); + boolList slavePCMask(saW.size(), false); + + forAll (partiallyCoveredMasterFaces, pfmI) + { + masterPCMask[partiallyCoveredMasterFaces[pfmI]] = true; + } + forAll (partiallyCoveredSlaveFaces, pfsI) + { + slavePCMask[partiallyCoveredSlaveFaces[pfsI]] = true; + } + // Rescaling the slave weights if (debug) { @@ -679,7 +746,7 @@ void GGIInterpolation::rescaleWeightingFactors() const { scalar slaveWeightSum = Foam::sum(saW[saWi]); - if (saW[saWi].size() > 0) + if (saW[saWi].size() > 0 && !slavePCMask[saWi]) { saW[saWi] = saW[saWi]/slaveWeightSum; @@ -696,7 +763,7 @@ void GGIInterpolation::rescaleWeightingFactors() const { scalar masterWeightSum = Foam::sum(maW[maWi]); - if (maW[maWi].size() > 0) + if (maW[maWi].size() > 0 && !masterPCMask[maWi]) { maW[maWi] = maW[maWi]/masterWeightSum; @@ -767,6 +834,119 @@ GGIInterpolation::findNonOverlappingFaces return tpatchFaceNonOverlapAddr; } + +template +void GGIInterpolation::calcPartiallyCoveredFaces +( + const scalarListList& patchWeights, + const scalar& nonOverlapFaceTol, + const bool isMaster +) const +{ + // Sanity checks first + if + ( + isMaster + && (partiallyCoveredMasterAddrPtr_ || masterFaceCoveredFractionsPtr_) + ) + { + FatalErrorIn + ( + "void GGIInterpolation::" + "calcPartiallyCoveredFaces() const" + ) << "Already calculated master partially covered faces" + << abort(FatalError); + } + + if + ( + !isMaster + && (partiallyCoveredSlaveAddrPtr_ || slaveFaceCoveredFractionsPtr_) + ) + { + FatalErrorIn + ( + "void GGIInterpolation::" + "calcPartiallyCoveredFaces() const" + ) << "Already calculated slave partially covered faces" + << abort(FatalError); + } + + // Temporary storage + DynamicList patchFacePartialOverlap(patchWeights.size()); + DynamicList uncoveredFaceFractions(patchWeights.size()); + + // Scan the list of patch weights and collect ones inbetween + // nonOverlapFaceTol and uncoveredFaceAreaTol_ + forAll (patchWeights, paWi) + { + const scalar sumWeightsFace = sum(patchWeights[paWi]); + + if + ( + sumWeightsFace > nonOverlapFaceTol + && sumWeightsFace <= uncoveredFaceAreaTol_() + ) + { + // This is considered partially overlapped face, store the index and + // the non-overlapping area (1 - sum of weights) + patchFacePartialOverlap.append(paWi); + uncoveredFaceFractions.append(1.0 - sumWeightsFace); + } + } + + // Transfer the storage + if (isMaster) + { + // Allocate master side + partiallyCoveredMasterAddrPtr_ = + new labelList(patchFacePartialOverlap.xfer()); + masterFaceCoveredFractionsPtr_ = + new scalarField(uncoveredFaceFractions.xfer()); + + if (debug) + { + InfoIn("GGIInterpolation::calcPartiallyCoveredFaces") + << " : Found " << partiallyCoveredMasterAddrPtr_->size() + << " partially overlapping faces for master GGI patch" << endl; + + if (partiallyCoveredMasterAddrPtr_->size()) + { + Info<< "Max coverage: " + << max(*masterFaceCoveredFractionsPtr_) + << ", min coverage: " + << min(*masterFaceCoveredFractionsPtr_) + << endl; + } + } + } + else + { + // Allocate slave side + partiallyCoveredSlaveAddrPtr_ = + new labelList(patchFacePartialOverlap.xfer()); + slaveFaceCoveredFractionsPtr_ = + new scalarField(uncoveredFaceFractions.xfer()); + + if (debug) + { + InfoIn("GGIInterpolation::calcPartiallyCoveredFaces") + << " : Found " << partiallyCoveredSlaveAddrPtr_->size() + << " partially overlapping faces for slave GGI patch" << endl; + + if (partiallyCoveredSlaveAddrPtr_->size()) + { + Info<< "Max coverage: " + << max(*slaveFaceCoveredFractionsPtr_) + << ", min coverage: " + << min(*slaveFaceCoveredFractionsPtr_) + << endl; + } + } + } +} + + template void GGIInterpolation:: calcMasterPointAddressing() const diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C index c9f6bc334..630572d89 100644 --- a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C +++ b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C @@ -324,7 +324,7 @@ void Foam::ggiPolyPatch::calcLocalParallel() const { FatalErrorIn("void ggiPolyPatch::calcLocalParallel() const") << "Patch size is greater than zone size for GGI patch " - << name() << ". This is not allowerd: " + << name() << ". This is not allowed: " << "the face zone must contain all patch faces and be " << "global in parallel runs" << abort(FatalError);