Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/compressible/rhoSimpleFoam/Allwclean

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wclean
+wclean rhoPorousMRFSimpleFoam
+wclean rhoSimplecFoam
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/compressible/rhoSimpleFoam/Allwmake

@@ -4,5 +4,6 @@ set -x
 
 wmake
 wmake rhoPorousMRFSimpleFoam
+wmake rhoSimplecFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/compressible/rhoSimpleFoam/UEqn.H

@@ -3,6 +3,7 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
+      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevRhoReff(U)
     );
 

applications/solvers/compressible/rhoSimpleFoam/hEqn.H

@@ -1,11 +1,13 @@
 {
+    volScalarField K("K", 0.5*magSqr(U));
+
     fvScalarMatrix hEqn
     (
         fvm::div(phi, h)
       - fvm::Sp(fvc::div(phi), h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-      - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
+        fvc::div(phi)*K - fvc::div(phi, K)
     );
 
     hEqn.relax();

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -28,7 +28,7 @@ if (simple.transonic())
         );
 
         // Relax the pressure equation to ensure diagonal-dominance
-        pEqn.relax(mesh.equationRelaxationFactor("pEqn"));
+        pEqn.relax();
 
         pEqn.setReference(pRefCell, pRefValue);
 

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options

@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I../rhoSimpleFoam \
+    -I.. \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H

@@ -29,21 +29,17 @@ if (simple.transonic())
         (
             "phic",
             fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf()
-          + phid*(fvc::interpolate(p) - fvc::interpolate(p, "UD"))
         );
 
-        //refCast<mixedFvPatchScalarField>(p.boundaryField()[1]).refValue()
-        //    = p.boundaryField()[1];
-
         fvScalarMatrix pEqn
         (
             fvm::div(phid, p)
           + fvc::div(phic)
-          - fvm::Sp(fvc::div(phid), p)
-          + fvc::div(phid)*p
           - fvm::laplacian(rho/AtU, p)
         );
-        //pEqn.relax();
+
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
 
         pEqn.setReference(pRefCell, pRefValue);
 
@@ -71,7 +67,6 @@ else
         fvScalarMatrix pEqn
         (
             fvc::div(phiHbyA)
-        //- fvm::laplacian(rho/AU, p)
           - fvm::laplacian(rho/AtU, p)
         );
 

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/lagrangian/sprayFoam/Allwmake

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake
+
+wmake sprayEngineFoam
+
+# ----------------------------------------------------------------- end-of-file

applications/test/tetTetOverlap/Make/files

@@ -0,0 +1,3 @@
+Test-tetTetOverlap.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-tetTetOverlap

applications/test/tetTetOverlap/Make/options

@@ -0,0 +1,5 @@
+EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools

applications/test/tetTetOverlap/Test-tetTetOverlap.C

@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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[])
+{
+    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(0.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;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/cvMesh/Allwmake

@@ -4,7 +4,7 @@ set -x
 
 wmake libso conformalVoronoiMesh
 wmake
-wmake cvMeshSurfaceSimplify
 wmake cvMeshBackgroundMesh
+(cd cvMeshSurfaceSimplify && ./Allwmake)
 
 # ----------------------------------------------------------------- end-of-file

applications/utilities/mesh/generation/cvMesh/Make/options

@@ -2,11 +2,15 @@ EXE_DEBUG = -DFULLDEBUG -g -O0
 EXE_FROUNDING_MATH = -frounding-math
 EXE_NDEBUG = -DNDEBUG
 
+CGAL_EXACT =
+CGAL_INEXACT = -DCGAL_INEXACT
+
 include $(GENERAL_RULES)/CGAL
 
 EXE_INC = \
     ${EXE_FROUNDING_MATH} \
     ${EXE_NDEBUG} \
+    ${CGAL_INEXACT} \
     ${CGAL_INC} \
     -IconformalVoronoiMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options

@@ -2,12 +2,16 @@ EXE_DEBUG = -DFULLDEBUG -g -O0
 EXE_FROUNDING_MATH = -frounding-math
 EXE_NDEBUG = -DNDEBUG
 
+CGAL_EXACT =
+CGAL_INEXACT = -DCGAL_INEXACT
+
 include $(GENERAL_RULES)/CGAL
 FFLAGS = -DCGAL_FILES='"${CGAL_ARCH_PATH}/share/files"'
 
 EXE_INC = \
     ${EXE_FROUNDING_MATH} \
     ${EXE_NDEBUG} \
+    ${CGAL_INEXACT} \
     ${CGAL_INC} \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C

@@ -249,6 +249,7 @@ Foam::conformationSurfaces::conformationSurfaces
 
     Info<< endl
         << "Testing for locationInMesh " << locationInMesh_ << endl;
+
     forAll(surfaces_, s)
     {
         const searchableSurface& surface(allGeometry_[surfaces_[s]]);
@@ -267,7 +268,6 @@ Foam::conformationSurfaces::conformationSurfaces
 
             referenceVolumeTypes_[s] = vTypes[0];
 
-
             Info<< "    is "
                 << searchableSurface::volumeTypeNames[referenceVolumeTypes_[s]]
                 << " surface " << surface.name()

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H

@@ -303,6 +303,9 @@ public:
             //- Find which patch is closest to the point
             label findPatch(const point& pt) const;
 
+            //- Is the surface a baffle
+            inline bool isBaffle(const label index) const;
+
 
         // Write
 

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfacesI.H

@@ -62,4 +62,10 @@ const Foam::treeBoundBox& Foam::conformationSurfaces::globalBounds() const
 }
 
 
+bool Foam::conformationSurfaces::isBaffle(const label index) const
+{
+    return baffleSurfaces_[index];
+}
+
+
 // ************************************************************************* //

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C

@@ -61,7 +61,7 @@ Foam::cvControls::cvControls
 
     specialiseFeaturePoints_ = Switch
     (
-        surfDict.lookupOrDefault<Switch>("specialiseFeaturePoints", false)
+        surfDict.lookup("specialiseFeaturePoints")
     );
 
     surfaceSearchDistanceCoeff_ = readScalar

applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/Allwmake

@@ -0,0 +1,10 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+if [ -d "${FASTDUALOCTREE_SRC_PATH}" ]
+then
+    wmake
+fi
+
+# ----------------------------------------------------------------- end-of-file

applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/Make/files

@@ -1,6 +1,22 @@
-cvMeshSurfaceSimplify.C
+/*
-
+cvMeshSurfaceSimplify_non_octree.C
 MarchingCubes/MarchingCubes.cpp
 MarchingCubes/ply.c
+*/
+
+/*
+MarchingCubes = fastdualoctree_sgp
+
+$(MarchingCubes)/data_access.cpp
+$(MarchingCubes)/fparser.cpp
+$(MarchingCubes)/fpoptimizer.cpp
+$(MarchingCubes)/MarchingCubes.cpp
+$(MarchingCubes)/mc_draw.cpp
+$(MarchingCubes)/morton.cpp
+$(MarchingCubes)/opt_octree.cpp
+$(MarchingCubes)/hash_octree.cpp
+*/
+
+cvMeshSurfaceSimplify.C
 
 EXE = $(FOAM_APPBIN)/cvMeshSurfaceSimplify

applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/Make/options

@@ -1,7 +1,12 @@
+MarchingCubes = fastdualoctree_sgp
+
 include $(GENERAL_RULES)/CGAL
 
 EXE_INC = \
-    -IMarchingCubes \
+    -DUNIX \
+    -Wno-old-style-cast \
+    /* -IMarchingCubes */ \
+    -I$(FASTDUALOCTREE_SRC_PATH) \
     -I../conformalVoronoiMesh/lnInclude \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
@@ -9,6 +14,8 @@ EXE_INC = \
 
 EXE_LIBS = \
     $(CGAL_LIBS) \
+    -L$(FASTDUALOCTREE_SRC_PATH) -lperf_main \
+    -lGL \
     -lconformalVoronoiMesh \
     -ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lscotchDecomp \
     -ledgeMesh \

applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/MarchingCubes/MarchingCubes.h

@@ -1,343 +0,0 @@
-/**
- * @file    MarchingCubes.h
- * @author  Thomas Lewiner <thomas.lewiner@polytechnique.org>
- * @author  Math Dept, PUC-Rio
- * @version 0.2
- * @date    12/08/2002
- *
- * @brief   MarchingCubes Algorithm
- */
-//________________________________________________
-
-
-#ifndef _MARCHINGCUBES_H_
-#define _MARCHINGCUBES_H_
-
-#if !defined(WIN32) || defined(__CYGWIN__)
-#pragma interface
-#endif // WIN32
-
-
-//_____________________________________________________________________________
-// types
-/** unsigned char alias */
-typedef unsigned char uchar ;
-/** signed char alias */
-typedef   signed char schar ;
-/** isovalue alias */
-typedef        float real  ;
-
-//-----------------------------------------------------------------------------
-// Vertex structure
-/** \struct Vertex "MarchingCubes.h" MarchingCubes
- * Position and normal of a vertex
- * \brief vertex structure
- * \param x X coordinate
- * \param y Y coordinate
- * \param z Z coordinate
- * \param nx X component of the normal
- * \param ny Y component of the normal
- * \param nz Z component of the normal
- */
-typedef struct
-{
-  real  x,  y,  z ;  /**< Vertex coordinates */
-  real nx, ny, nz ;  /**< Vertex normal */
-} Vertex ;
-
-//-----------------------------------------------------------------------------
-// Triangle structure
-/** \struct Triangle "MarchingCubes.h" MarchingCubes
- * Indices of the oriented triange vertices
- * \brief triangle structure
- * \param v1 First vertex index
- * \param v2 Second vertex index
- * \param v3 Third vertex index
- */
-typedef struct
-{
-  int v1,v2,v3 ;  /**< Triangle vertices */
-} Triangle ;
-//_____________________________________________________________________________
-
-
-
-//_____________________________________________________________________________
-/** Marching Cubes algorithm wrapper */
-/** \class MarchingCubes
-  * \brief Marching Cubes algorithm.
-  */
-class MarchingCubes
-//-----------------------------------------------------------------------------
-{
-// Constructors
-public :
-  /**
-   * Main and default constructor
-   * \brief constructor
-   * \param size_x width  of the grid
-   * \param size_y depth  of the grid
-   * \param size_z height of the grid
-   */
-  MarchingCubes ( const int size_x = -1, const int size_y = -1, const int size_z = -1 ) ;
-  /** Destructor */
-  ~MarchingCubes() ;
-
-//-----------------------------------------------------------------------------
-// Accessors
-public :
-  /** accesses the number of vertices of the generated mesh */
-  inline const int nverts() const { return _nverts ; }
-  /** accesses the number of triangles of the generated mesh */
-  inline const int ntrigs() const { return _ntrigs ; }
-  /** accesses a specific vertex of the generated mesh */
-  inline Vertex   * vert( const int i ) const { if( i < 0  || i >= _nverts ) return ( Vertex *)NULL ; return _vertices  + i ; }
-  /** accesses a specific triangle of the generated mesh */
-  inline Triangle * trig( const int i ) const { if( i < 0  || i >= _ntrigs ) return (Triangle*)NULL ; return _triangles + i ; }
-
-  /** accesses the vertex buffer of the generated mesh */
-  inline Vertex   *vertices () { return _vertices  ; }
-  /** accesses the triangle buffer of the generated mesh */
-  inline Triangle *triangles() { return _triangles ; }
-
-  /**  accesses the width  of the grid */
-  inline const int size_x() const { return _size_x ; }
-  /**  accesses the depth  of the grid */
-  inline const int size_y() const { return _size_y ; }
-  /**  accesses the height of the grid */
-  inline const int size_z() const { return _size_z ; }
-
-  /**
-   * changes the size of the grid
-   * \param size_x width  of the grid
-   * \param size_y depth  of the grid
-   * \param size_z height of the grid
-   */
-  inline void set_resolution( const int size_x, const int size_y, const int size_z ) { _size_x = size_x ;  _size_y = size_y ;  _size_z = size_z ; }
-  /**
-   * selects wether the algorithm will use the enhanced topologically controlled lookup table or the original MarchingCubes
-   * \param originalMC true for the original Marching Cubes
-   */
-  inline void set_method    ( const bool originalMC = false ) { _originalMC = originalMC ; }
-  /**
-   * selects to use data from another class
-   * \param data is the pointer to the external data, allocated as a size_x*size_y*size_z vector running in x first
-   */
-  inline void set_ext_data  ( real *data )
-  { if( !_ext_data ) delete [] _data ;  _ext_data = data != NULL ;  if( _ext_data ) _data = data ; }
-  /**
-   * selects to allocate data
-   */
-  inline void set_int_data  () { _ext_data = false ;  _data = NULL ; }
-
-  // Data access
-  /**
-   * accesses a specific cube of the grid
-   * \param i abscisse of the cube
-   * \param j ordinate of the cube
-   * \param k height of the cube
-   */
-  inline const real get_data  ( const int i, const int j, const int k ) const { return _data[ i + j*_size_x + k*_size_x*_size_y] ; }
-  /**
-   * sets a specific cube of the grid
-   * \param val new value for the cube
-   * \param i abscisse of the cube
-   * \param j ordinate of the cube
-   * \param k height of the cube
-   */
-  inline void  set_data  ( const real val, const int i, const int j, const int k ) { _data[ i + j*_size_x + k*_size_x*_size_y] = val ; }
-
-  // Data initialization
-  /** inits temporary structures (must set sizes before call) : the grid and the vertex index per cube */
-  void init_temps () ;
-  /** inits all structures (must set sizes before call) : the temporary structures and the mesh buffers */
-  void init_all   () ;
-  /** clears temporary structures : the grid and the main */
-  void clean_temps() ;
-  /** clears all structures : the temporary structures and the mesh buffers */
-  void clean_all  () ;
-
-
-//-----------------------------------------------------------------------------
-// Exportation
-public :
-  /**
-   * PLY exportation of the generated mesh
-   * \param fn  name of the PLY file to create
-   * \param bin if true, the PLY will be written in binary mode
-   */
-  void writePLY( const char *fn, bool bin = false ) ;
-
-  /**
-   * PLY importation of a mesh
-   * \param fn  name of the PLY file to read from
-   */
-  void readPLY( const char *fn ) ;
-
-  /**
-   * VRML / Open Inventor exportation of the generated mesh
-   * \param fn  name of the IV file to create
-   */
-  void writeIV ( const char *fn ) ;
-
-  /**
-   * ISO exportation of the input grid
-   * \param fn  name of the ISO file to create
-   */
-  void writeISO( const char *fn ) ;
-
-
-//-----------------------------------------------------------------------------
-// Algorithm
-public :
-  /**
-   * Main algorithm : must be called after init_all
-   * \param iso isovalue
-   */
-  void run( real iso = (real)0.0 ) ;
-
-protected :
-  /** tesselates one cube */
-  void process_cube ()             ;
-  /** tests if the components of the tesselation of the cube should be connected by the interior of an ambiguous face */
-  bool test_face    ( schar face ) ;
-  /** tests if the components of the tesselation of the cube should be connected through the interior of the cube */
-  bool test_interior( schar s )    ;
-
-
-//-----------------------------------------------------------------------------
-// Operations
-protected :
-  /**
-   * computes almost all the vertices of the mesh by interpolation along the cubes edges
-   * \param iso isovalue
-   */
-  void compute_intersection_points( real iso ) ;
-
-  /**
-   * routine to add a triangle to the mesh
-   * \param trig the code for the triangle as a sequence of edges index
-   * \param n    the number of triangles to produce
-   * \param v12  the index of the interior vertex to use, if necessary
-   */
-  void add_triangle ( const char* trig, char n, int v12 = -1 ) ;
-
-  /** tests and eventually doubles the vertex buffer capacity for a new vertex insertion */
-  void test_vertex_addition() ;
-  /** adds a vertex on the current horizontal edge */
-  int add_x_vertex() ;
-  /** adds a vertex on the current longitudinal edge */
-  int add_y_vertex() ;
-  /** adds a vertex on the current vertical edge */
-  int add_z_vertex() ;
-  /** adds a vertex inside the current cube */
-  int add_c_vertex() ;
-
-  /**
-   * interpolates the horizontal gradient of the implicit function at the lower vertex of the specified cube
-   * \param i abscisse of the cube
-   * \param j ordinate of the cube
-   * \param k height of the cube
-   */
-  real get_x_grad( const int i, const int j, const int k ) const ;
-  /**
-   * interpolates the longitudinal gradient of the implicit function at the lower vertex of the specified cube
-   * \param i abscisse of the cube
-   * \param j ordinate of the cube
-   * \param k height of the cube
-   */
-  real get_y_grad( const int i, const int j, const int k ) const ;
-  /**
-   * interpolates the vertical gradient of the implicit function at the lower vertex of the specified cube
-   * \param i abscisse of the cube
-   * \param j ordinate of the cube
-   * \param k height of the cube
-   */
-  real get_z_grad( const int i, const int j, const int k ) const ;
-
-  /**
-   * accesses the pre-computed vertex index on the lower horizontal edge of a specific cube
-   * \param i abscisse of the cube
-   * \param j ordinate of the cube
-   * \param k height of the cube
-   */
-  inline int   get_x_vert( const int i, const int j, const int k ) const { return _x_verts[ i + j*_size_x + k*_size_x*_size_y] ; }
-  /**
-   * accesses the pre-computed vertex index on the lower longitudinal edge of a specific cube
-   * \param i abscisse of the cube
-   * \param j ordinate of the cube
-   * \param k height of the cube
-   */
-  inline int   get_y_vert( const int i, const int j, const int k ) const { return _y_verts[ i + j*_size_x + k*_size_x*_size_y] ; }
-  /**
-   * accesses the pre-computed vertex index on the lower vertical edge of a specific cube
-   * \param i abscisse of the cube
-   * \param j ordinate of the cube
-   * \param k height of the cube
-   */
-  inline int   get_z_vert( const int i, const int j, const int k ) const { return _z_verts[ i + j*_size_x + k*_size_x*_size_y] ; }
-
-  /**
-   * sets the pre-computed vertex index on the lower horizontal edge of a specific cube
-   * \param val the index of the new vertex
-   * \param i abscisse of the cube
-   * \param j ordinate of the cube
-   * \param k height of the cube
-   */
-  inline void  set_x_vert( const int val, const int i, const int j, const int k ) { _x_verts[ i + j*_size_x + k*_size_x*_size_y] = val ; }
-  /**
-   * sets the pre-computed vertex index on the lower longitudinal edge of a specific cube
-   * \param val the index of the new vertex
-   * \param i abscisse of the cube
-   * \param j ordinate of the cube
-   * \param k height of the cube
-   */
-  inline void  set_y_vert( const int val, const int i, const int j, const int k ) { _y_verts[ i + j*_size_x + k*_size_x*_size_y] = val ; }
-  /**
-   * sets the pre-computed vertex index on the lower vertical edge of a specific cube
-   * \param val the index of the new vertex
-   * \param i abscisse of the cube
-   * \param j ordinate of the cube
-   * \param k height of the cube
-   */
-  inline void  set_z_vert( const int val, const int i, const int j, const int k ) { _z_verts[ i + j*_size_x + k*_size_x*_size_y] = val ; }
-
-  /** prints cube for debug */
-  void print_cube() ;
-
-//-----------------------------------------------------------------------------
-// Elements
-protected :
-  bool      _originalMC ;   /**< selects wether the algorithm will use the enhanced topologically controlled lookup table or the original MarchingCubes */
-  bool      _ext_data   ;   /**< selects wether to allocate data or use data from another class */
-
-  int       _size_x     ;  /**< width  of the grid */
-  int       _size_y     ;  /**< depth  of the grid */
-  int       _size_z     ;  /**< height of the grid */
-  real     *_data       ;  /**< implicit function values sampled on the grid */
-
-  int      *_x_verts    ;  /**< pre-computed vertex indices on the lower horizontal   edge of each cube */
-  int      *_y_verts    ;  /**< pre-computed vertex indices on the lower longitudinal edge of each cube */
-  int      *_z_verts    ;  /**< pre-computed vertex indices on the lower vertical     edge of each cube */
-
-  int       _nverts     ;  /**< number of allocated vertices  in the vertex   buffer */
-  int       _ntrigs     ;  /**< number of allocated triangles in the triangle buffer */
-  int       _Nverts     ;  /**< size of the vertex   buffer */
-  int       _Ntrigs     ;  /**< size of the triangle buffer */
-  Vertex   *_vertices   ;  /**< vertex   buffer */
-  Triangle *_triangles  ;  /**< triangle buffer */
-
-  int       _i          ;  /**< abscisse of the active cube */
-  int       _j          ;  /**< height of the active cube */
-  int       _k          ;  /**< ordinate of the active cube */
-
-  real      _cube[8]    ;  /**< values of the implicit function on the active cube */
-  uchar     _lut_entry  ;  /**< cube sign representation in [0..255] */
-  uchar     _case       ;  /**< case of the active cube in [0..15] */
-  uchar     _config     ;  /**< configuration of the active cube */
-  uchar     _subconfig  ;  /**< subconfiguration of the active cube */
-};
-//_____________________________________________________________________________
-
-
-#endif // _MARCHINGCUBES_H_

applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/MarchingCubes/ply.h

@@ -1,233 +0,0 @@
-/*
-
-Header for PLY polygon files.
-
-- Greg Turk
-
-A PLY file contains a single polygonal _object_.
-
-An object is composed of lists of _elements_.  Typical elements are
-vertices, faces, edges and materials.
-
-Each type of element for a given object has one or more _properties_
-associated with the element type.  For instance, a vertex element may
-have as properties three floating-point values x,y,z and three unsigned
-chars for red, green and blue.
-
------------------------------------------------------------------------
-
-Copyright (c) 1998 Georgia Institute of Technology.  All rights reserved.
-
-Permission to use, copy, modify and distribute this software and its
-documentation for any purpose is hereby granted without fee, provided
-that the above copyright notice and this permission notice appear in
-all copies of this software and that you do not sell the software.
-
-THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
-EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
-WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
-
-*/
-
-#ifndef __PLY_H__
-#define __PLY_H__
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-#include <stdio.h>
-#include <stddef.h>
-
-#define PLY_ASCII      1        /* ascii PLY file */
-#define PLY_BINARY_BE  2        /* binary PLY file, big endian */
-#define PLY_BINARY_LE  3        /* binary PLY file, little endian */
-
-#define PLY_OKAY    0           /* ply routine worked okay */
-#define PLY_ERROR  -1           /* error in ply routine */
-
-/* scalar data types supported by PLY format */
-
-#define StartType  0
-#define Int8       1
-#define Int16      2
-#define Int32      3
-#define Uint8      4
-#define Uint16     5
-#define Uint32     6
-#define Float32    7
-#define Float64    8
-#define EndType    9
-
-#define  PLY_SCALAR  0
-#define  PLY_LIST    1
-#define  PLY_STRING  2
-
-
-typedef struct PlyProperty {    /* description of a property */
-
-char *name;                   /* property name */
-int external_type;            /* file's data type */
-int internal_type;            /* program's data type */
-int offset;                   /* offset bytes of prop in a struct */
-
-int is_list;                  /* 0 = scalar, 1 = list, 2 = char string */
-int count_external;           /* file's count type */
-int count_internal;           /* program's count type */
-int count_offset;             /* offset byte for list count */
-
-} PlyProperty;
-
-typedef struct PlyElement {     /* description of an element */
-char *name;                   /* element name */
-int num;                      /* number of elements in this object */
-int size;                     /* size of element (bytes) or -1 if variable */
-int nprops;                   /* number of properties for this element */
-PlyProperty **props;          /* list of properties in the file */
-char *store_prop;             /* flags: property wanted by user? */
-int other_offset;             /* offset to un-asked-for props, or -1 if none*/
-int other_size;               /* size of other_props structure */
-} PlyElement;
-
-typedef struct PlyOtherProp {   /* describes other properties in an element */
-char *name;                   /* element name */
-int size;                     /* size of other_props */
-int nprops;                   /* number of properties in other_props */
-PlyProperty **props;          /* list of properties in other_props */
-} PlyOtherProp;
-
-typedef struct OtherData { /* for storing other_props for an other element */
-void *other_props;
-} OtherData;
-
-typedef struct OtherElem {     /* data for one "other" element */
-char *elem_name;             /* names of other elements */
-int elem_count;              /* count of instances of each element */
-OtherData **other_data;      /* actual property data for the elements */
-PlyOtherProp *other_props;   /* description of the property data */
-} OtherElem;
-
-typedef struct PlyOtherElems {  /* "other" elements, not interpreted by user */
-int num_elems;                /* number of other elements */
-OtherElem *other_list;        /* list of data for other elements */
-} PlyOtherElems;
-
-#define AVERAGE_RULE  1
-#define MAJORITY_RULE 2
-#define MINIMUM_RULE  3
-#define MAXIMUM_RULE  4
-#define SAME_RULE     5
-#define RANDOM_RULE   6
-
-typedef struct PlyPropRules {   /* rules for combining "other" properties */
-PlyElement *elem;      /* element whose rules we are making */
-int *rule_list;        /* types of rules (AVERAGE_PLY, MAJORITY_PLY, etc.) */
-int nprops;            /* number of properties we're combining so far */
-int max_props;         /* maximum number of properties we have room for now */
-void **props;          /* list of properties we're combining */
-float *weights;        /* list of weights of the properties */
-} PlyPropRules;
-
-typedef struct PlyRuleList {
-char *name;                  /* name of the rule */
-char *element;               /* name of element that rule applies to */
-char *property;              /* name of property that rule applies to */
-struct PlyRuleList *next;    /* pointer for linked list of rules */
-} PlyRuleList;
-
-typedef struct PlyFile {        /* description of PLY file */
-FILE *fp;                     /* file pointer */
-int file_type;                /* ascii or binary */
-float version;                /* version number of file */
-int num_elem_types;           /* number of element types of object */
-PlyElement **elems;           /* list of elements */
-int num_comments;             /* number of comments */
-char **comments;              /* list of comments */
-int num_obj_info;             /* number of items of object information */
-char **obj_info;              /* list of object info items */
-PlyElement *which_elem;       /* element we're currently reading or writing */
-PlyOtherElems *other_elems;   /* "other" elements from a PLY file */
-PlyPropRules *current_rules;  /* current propagation rules */
-PlyRuleList *rule_list;       /* rule list from user */
-} PlyFile;
-
-/* memory allocation */
-/*
-extern char *my_alloc();
-*/
-#define myalloc(mem_size) my_alloc((mem_size), __LINE__, __FILE__)
-
-
-/* old routines */
-
-#if 0
-extern PlyFile *ply_write(FILE *, int, char **, int);
-extern PlyFile *ply_read(FILE *, int *, char ***);
-extern PlyFile *ply_open_for_reading( char *, int *, char ***, int *, float *);
-extern void ply_close(PlyFile *);
-extern PlyOtherProp *ply_get_other_properties(PlyFile *, char *, int);
-#endif
-
-extern void     ply_describe_property( PlyFile * , char * , PlyProperty * );
-extern void     ply_get_property( PlyFile * , char * , PlyProperty * );
-extern void     ply_get_element( PlyFile * , void * );
-
-
-/*** delcaration of routines ***/
-
-PlyOtherElems  *get_other_element_ply( PlyFile * );
-
-PlyFile        *read_ply( FILE * );
-PlyFile        *write_ply( FILE * , int, char ** , int );
-extern PlyFile *open_for_writing_ply( char * , int, char ** , int );
-void            close_ply( PlyFile * );
-void            free_ply( PlyFile * );
-
-void            get_info_ply( PlyFile * , float * , int * );
-void            free_other_elements_ply( PlyOtherElems * );
-
-void            append_comment_ply( PlyFile * , char * );
-void            append_obj_info_ply( PlyFile * , char * );
-void            copy_comments_ply( PlyFile * , PlyFile * );
-void            copy_obj_info_ply( PlyFile * , PlyFile * );
-char*          *get_comments_ply( PlyFile * , int * );
-char*          *get_obj_info_ply( PlyFile * , int * );
-
-char*          *get_element_list_ply( PlyFile * , int * );
-int             setup_property_ply( PlyFile * , PlyProperty * );
-void            get_element_ply( PlyFile * , void * );
-char           *setup_element_read_ply( PlyFile * , int, int * );
-PlyOtherProp   *get_other_properties_ply( PlyFile * , int );
-
-void            element_count_ply( PlyFile * , char * , int );
-void            describe_element_ply( PlyFile * , char * , int );
-void            describe_property_ply( PlyFile * , PlyProperty * );
-void            describe_other_properties_ply( PlyFile * , PlyOtherProp * , int );
-void            describe_other_elements_ply( PlyFile * , PlyOtherElems * );
-void            get_element_setup_ply( PlyFile * , char * , int, PlyProperty * );
-PlyProperty*   *get_element_description_ply( PlyFile * , char * , int * , int * );
-void            element_layout_ply( PlyFile * , char * , int, int, PlyProperty * );
-
-void            header_complete_ply( PlyFile * );
-void            put_element_setup_ply( PlyFile * , char * );
-void            put_element_ply( PlyFile * , void * );
-void            put_other_elements_ply( PlyFile * );
-
-PlyPropRules   *init_rule_ply( PlyFile * , char * );
-void            modify_rule_ply( PlyPropRules * , char * , int );
-void            start_props_ply( PlyFile * , PlyPropRules * );
-void            weight_props_ply( PlyFile * , float, void * );
-void           *get_new_props_ply( PlyFile * );
-void            set_prop_rules_ply( PlyFile * , PlyRuleList * );
-PlyRuleList    *append_prop_rule( PlyRuleList * , char * , char * );
-int             matches_rule_name( char * );
-
-int             equal_strings( char * , char * );
-char           *recreate_command_line( int, char *argv[] );
-
-
-#ifdef __cplusplus
-}
-#endif
-#endif /* !__PLY_H__ */
-

applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/cvMeshSurfaceSimplify.C

@@ -28,6 +28,7 @@ Description
     Simplifies surfaces by resampling.
 
     Uses Thomas Lewiner's topology preserving MarchingCubes.
+    (http://zeus.mat.puc-rio.br/tomlew/tomlew_uk.php)
 
 \*---------------------------------------------------------------------------*/
 
@@ -37,13 +38,328 @@ Description
 #include "conformationSurfaces.H"
 #include "triSurfaceMesh.H"
 
-#include "MarchingCubes.h"
+#include "opt_octree.h"
-
+#include "cube.h"
 
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+class pointConversion
+{
+    const vector scale_;
+
+    const vector offset_;
+
+public:
+
+    //- Construct from components
+    pointConversion
+    (
+        const vector scale,
+        const vector offset
+    )
+    :
+        scale_(scale),
+        offset_(offset)
+    {}
+
+    inline Point toLocal(const Foam::point& pt) const
+    {
+        Foam::point p = cmptMultiply(scale_, (pt + offset_));
+        return Point(p.x(), p.y(), p.z());
+    }
+
+    inline Foam::point toGlobal(const Point& pt) const
+    {
+        point p(pt.x(), pt.y(), pt.z());
+        return cmptDivide(p, scale_) - offset_;
+    }
+};
+
+
+
+
+// For use in Fast-Dual Octree from Thomas Lewiner
+class distanceCalc
+:
+    public ::data_access
+{
+
+    const Level min_level_;
+
+    const conformationSurfaces& geometryToConformTo_;
+
+    const pointConversion& converter_;
+
+
+    // Private Member Functions
+
+    scalar signedDistance(const Foam::point& pt) const
+    {
+        const searchableSurfaces& geometry = geometryToConformTo_.geometry();
+        const labelList& surfaces = geometryToConformTo_.surfaces();
+
+        static labelList nearestSurfaces;
+        static scalarField distance;
+
+        static pointField samples(1);
+        samples[0] = pt;
+
+        searchableSurfacesQueries::signedDistance
+        (
+            geometry,
+            surfaces,
+            samples,
+            scalarField(1, GREAT),
+            searchableSurface::OUTSIDE,
+            nearestSurfaces,
+            distance
+        );
+
+        return distance[0];
+    }
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        distanceCalc
+        (
+            Level max_level_,
+            real iso_val_,
+            Level min_level,
+            const conformationSurfaces& geometryToConformTo,
+            const pointConversion& converter
+        )
+        :
+            data_access(max_level_,iso_val_),
+            min_level_(min_level),
+            geometryToConformTo_(geometryToConformTo),
+            converter_(converter)
+        {}
+
+
+    //- Destructor
+    virtual ~distanceCalc()
+    {}
+
+
+    // Member Functions
+
+        //- test function
+        virtual bool need_refine( const Cube &c )
+        {
+            int l = c.lv() ;
+
+            if( l >= _max_level ) return false;
+            if( l < min_level_ ) return true;
+
+            treeBoundBox bb
+            (
+                converter_.toGlobal
+                (
+                    Point
+                    (
+                        c.xmin(),
+                        c.ymin(),
+                        c.zmin()
+                    )
+                ),
+                converter_.toGlobal
+                (
+                    Point
+                    (
+                        c.xmax(),
+                        c.ymax(),
+                        c.zmax()
+                    )
+                )
+            );
+
+            const searchableSurfaces& geometry =
+                geometryToConformTo_.geometry();
+            const labelList& surfaces =
+                geometryToConformTo_.surfaces();
+
+
+            //- Uniform refinement around surface
+            {
+                forAll(surfaces, i)
+                {
+                    if (geometry[surfaces[i]].overlaps(bb))
+                    {
+                        return true;
+                    }
+                }
+                return false;
+            }
+
+
+            ////- Surface intersects bb (but not using intersection test)
+            //scalar ccDist = signedDistance(bb.midpoint());
+            //scalar ccVal = ccDist - _iso_val;
+            //if (mag(ccVal) < SMALL)
+            //{
+            //    return true;
+            //}
+            //const pointField points(bb.points());
+            //forAll(points, pointI)
+            //{
+            //    scalar pointVal = signedDistance(points[pointI]) - _iso_val;
+            //    if (ccVal*pointVal < 0)
+            //    {
+            //        return true;
+            //    }
+            //}
+            //return false;
+
+
+            ////- Refinement based on intersection with multiple planes.
+            ////  Does not work well - too high a ratio between
+            ////  neighbouring cubes.
+            //const pointField points(bb.points());
+            //const edgeList& edges = treeBoundBox::edges;
+            //pointField start(edges.size());
+            //pointField end(edges.size());
+            //forAll(edges, i)
+            //{
+            //    start[i] = points[edges[i][0]];
+            //    end[i] = points[edges[i][1]];
+            //}
+            //Foam::List<Foam::List<pointIndexHit> > hitInfo;
+            //labelListList hitSurfaces;
+            //searchableSurfacesQueries::findAllIntersections
+            //(
+            //    geometry,
+            //    surfaces,
+            //    start,
+            //    end,
+            //    hitSurfaces,
+            //    hitInfo
+            //);
+            //
+            //// Count number of intersections
+            //label nInt = 0;
+            //forAll(hitSurfaces, edgeI)
+            //{
+            //    nInt += hitSurfaces[edgeI].size();
+            //}
+            //
+            //if (nInt == 0)
+            //{
+            //    // No surface intersected. See if there is one inside
+            //    forAll(surfaces, i)
+            //    {
+            //        if (geometry[surfaces[i]].overlaps(bb))
+            //        {
+            //            return true;
+            //        }
+            //    }
+            //    return false;
+            //}
+            //
+            //// Check multiple surfaces
+            //label baseSurfI = -1;
+            //forAll(hitSurfaces, edgeI)
+            //{
+            //    const labelList& hSurfs = hitSurfaces[edgeI];
+            //    forAll(hSurfs, i)
+            //    {
+            //        if (baseSurfI == -1)
+            //        {
+            //            baseSurfI = hSurfs[i];
+            //        }
+            //        else if (baseSurfI != hSurfs[i])
+            //        {
+            //            // Multiple surfaces
+            //            return true;
+            //        }
+            //    }
+            //}
+            //
+            //// Get normals
+            //DynamicList<pointIndexHit> baseInfo(nInt);
+            //forAll(hitInfo, edgeI)
+            //{
+            //    const Foam::List<pointIndexHit>& hits = hitInfo[edgeI];
+            //    forAll(hits, i)
+            //    {
+            //        (void)hits[i].hitPoint();
+            //        baseInfo.append(hits[i]);
+            //    }
+            //}
+            //vectorField normals;
+            //geometry[surfaces[baseSurfI]].getNormal(baseInfo, normals);
+            //for (label i = 1; i < normals.size(); ++i)
+            //{
+            //    if ((normals[0] & normals[i]) < 0.9)
+            //    {
+            //        return true;
+            //    }
+            //}
+            //labelList regions;
+            //geometry[surfaces[baseSurfI]].getRegion(baseInfo, regions);
+            //for (label i = 1; i < regions.size(); ++i)
+            //{
+            //    if (regions[0] != regions[i])
+            //    {
+            //        return true;
+            //    }
+            //}
+            //return false;
+
+
+
+            //samples[0] = point(c.xmin(), c.ymin(), c.zmin());
+            //samples[1] = point(c.xmax(), c.ymin(), c.zmin());
+            //samples[2] = point(c.xmax(), c.ymax(), c.zmin());
+            //samples[3] = point(c.xmin(), c.ymax(), c.zmin());
+            //
+            //samples[4] = point(c.xmin(), c.ymin(), c.zmax());
+            //samples[5] = point(c.xmax(), c.ymin(), c.zmax());
+            //samples[6] = point(c.xmax(), c.ymax(), c.zmax());
+            //samples[7] = point(c.xmin(), c.ymax(), c.zmax());
+
+            //scalarField nearestDistSqr(8, GREAT);
+            //
+            //Foam::List<pointIndexHit> nearestInfo;
+            //surf_.findNearest(samples, nearestDistSqr, nearestInfo);
+            //vectorField normals;
+            //surf_.getNormal(nearestInfo, normals);
+            //
+            //for (label i = 1; i < normals.size(); ++i)
+            //{
+            //    if ((normals[0] & normals[i]) < 0.5)
+            //    {
+            //        return true;
+            //    }
+            //}
+            //return false;
+
+            //// Check if surface octree same level
+            //const labelList elems(surf_.tree().findBox(bb));
+            //
+            //if (elems.size() > 1)
+            //{
+            //    return true;
+            //}
+            //else
+            //{
+            //  return false;
+            //}
+        }
+
+        //- data function
+        virtual real value_at( const Cube &c )
+        {
+            return signedDistance(converter_.toGlobal(c)) - _iso_val;
+        }
+};
+
+
 // Main program:
 
 int main(int argc, char *argv[])
@@ -52,19 +368,16 @@ int main(int argc, char *argv[])
     (
         "Re-sample surfaces used in cvMesh operation"
     );
-    //argList::validArgs.append("inputFile");
-    argList::validArgs.append("(nx ny nz)");
     argList::validArgs.append("outputName");
 
     #include "setRootCase.H"
     #include "createTime.H"
     runTime.functionObjects().off();
 
-    const Vector<label> n(IStringStream(args.args()[1])());
+    const fileName exportName = args.args()[1];
-    const fileName exportName = args.args()[2];
 
     Info<< "Reading surfaces as specified in the cvMeshDict and"
-        << " writing re-sampled " << n << " to " << exportName
+        << " writing a re-sampled surface to " << exportName
         << nl << endl;
 
     cpuTime timer;
@@ -114,126 +427,100 @@ int main(int argc, char *argv[])
         << timer.cpuTimeIncrement() << " s." << nl << endl;
 
 
-
-    // Extend
-    treeBoundBox bb = geometryToConformTo.globalBounds();
-    {
-        const vector smallVec = 0.1*bb.span();
-        bb.min() -= smallVec;
-        bb.max() += smallVec;
-    }
-
-    Info<< "Meshing to bounding box " << bb << nl << endl;
-
-    const vector span(bb.span());
-    const vector d
-    (
-        span.x()/(n.x()-1),
-        span.y()/(n.y()-1),
-        span.z()/(n.z()-1)
-    );
-
-    MarchingCubes mc(span.x(), span.y(), span.z() ) ;
-    mc.set_resolution(n.x(), n.y(), n.z());
-    mc.init_all() ;
-
-
-    // Generate points
-    pointField points(mc.size_x()*mc.size_y()*mc.size_z());
-    label pointI = 0;
-
-    point pt;
-    for( int k = 0 ; k < mc.size_z() ; k++ )
-    {
-        pt.z() = bb.min().z() + k*d.z();
-        for( int j = 0 ; j < mc.size_y() ; j++ )
-        {
-            pt.y() = bb.min().y() + j*d.y();
-            for( int i = 0 ; i < mc.size_x() ; i++ )
-            {
-                pt.x() = bb.min().x() + i*d.x();
-                points[pointI++] = pt;
-            }
-        }
-    }
-
-    Info<< "Generated " << points.size() << " sampling points in = "
-        << timer.cpuTimeIncrement() << " s." << nl << endl;
-
-
-    // Compute data
     const searchableSurfaces& geometry = geometryToConformTo.geometry();
     const labelList& surfaces = geometryToConformTo.surfaces();
 
-    scalarField signedDist;
+
-    labelList nearestSurfaces;
+    const label minLevel = 2;
-    searchableSurfacesQueries::signedDistance
+
+    // The max cube size follows from the minLevel and the default cube size
+    // (1)
+    const scalar maxSize = 1.0 / (1 << minLevel);
+    const scalar halfMaxSize = 0.5*maxSize;
+
+
+    // Scale the geometry to fit within
+    // halfMaxSize .. 1-halfMaxSize
+
+    scalar wantedRange = 1.0-maxSize;
+
+    const treeBoundBox bb = geometryToConformTo.globalBounds();
+
+    const vector scale = cmptDivide
     (
-        geometry,
+        point(wantedRange, wantedRange, wantedRange),
-        surfaces,
+        bb.span()
-        points,
+    );
-        scalarField(points.size(), sqr(GREAT)),
+    const vector offset =
-        searchableSurface::OUTSIDE,     // for non-closed surfaces treat as
+        cmptDivide
-                                        // outside
+        (
-        nearestSurfaces,
+            point(halfMaxSize, halfMaxSize, halfMaxSize),
-        signedDist
+            scale
+        )
+       -bb.min();
+
+
+    const pointConversion converter(scale, offset);
+
+
+    // Marching cubes
+
+    OptOctree octree;
+
+    distanceCalc ref
+    (
+        8,          //maxLevel
+        0.0,        //distance
+        minLevel,   //minLevel
+        geometryToConformTo,
+        converter
     );
 
-    // Fill elements
+    octree.refine(&ref);
-    pointI = 0;
+    octree.set_impl(&ref);
-    for( int k = 0 ; k < mc.size_z() ; k++ )
-    {
-        for( int j = 0 ; j < mc.size_y() ; j++ )
-        {
-            for( int i = 0 ; i < mc.size_x() ; i++ )
-            {
-                mc.set_data(float(signedDist[pointI++]), i, j, k);
-            }
-        }
-    }
 
-    Info<< "Determined signed distance in = "
+    Info<< "Calculated octree in = "
         << timer.cpuTimeIncrement() << " s." << nl << endl;
 
+    MarchingCubes& mc = octree.mc();
 
-    mc.run() ;
+    mc.clean_all() ;
+    octree.build_isosurface(&ref) ;
 
-    Info<< "Constructed iso surface in = "
+    Info<< "Constructed iso surface of distance in = "
         << timer.cpuTimeIncrement() << " s." << nl << endl;
 
-
-    mc.clean_temps() ;
-
-
-
     // Write output file
     if (mc.ntrigs() > 0)
     {
         Triangle* triangles = mc.triangles();
-        List<labelledTri> tris(mc.ntrigs());
+        label nTris = mc.ntrigs();
-        forAll(tris, triI)
+        Foam::DynamicList<labelledTri> tris(mc.ntrigs());
+        for (label triI = 0; triI < nTris; ++triI)
         {
-            tris[triI] = labelledTri
+            const Triangle& t = triangles[triI];
-            (
+            if (t.v1 != t.v2 && t.v1 != t.v3 && t.v2 != t.v3)
-                triangles[triI].v1,
+            {
-                triangles[triI].v2,
+                tris.append
-                triangles[triI].v3,
+                (
-                0                       // region
+                    labelledTri
-            );
+                    (
+                        triangles[triI].v1,
+                        triangles[triI].v2,
+                        triangles[triI].v3,
+                        0                       // region
+                    )
+                );
+            }
         }
 
 
-        Vertex* vertices = mc.vertices();
+        Point* vertices = mc.vertices();
         pointField points(mc.nverts());
         forAll(points, pointI)
         {
-            Vertex& v = vertices[pointI];
+            const Point& v = vertices[pointI];
-            points[pointI] = point
+            points[pointI] = converter.toGlobal(v);
-            (
-                bb.min().x() + v.x*span.x()/mc.size_x(),
-                bb.min().y() + v.y*span.y()/mc.size_y(),
-                bb.min().z() + v.z*span.z()/mc.size_z()
-            );
         }
 
 
@@ -267,6 +554,7 @@ int main(int argc, char *argv[])
         }
 
         triSurface s(tris, patches, points, true);
+        tris.clearStorage();
 
         Info<< "Extracted triSurface in = "
             << timer.cpuTimeIncrement() << " s." << nl << endl;
@@ -274,7 +562,7 @@ int main(int argc, char *argv[])
 
         // Find out region on local surface of nearest point
         {
-            List<pointIndexHit> hitInfo;
+            Foam::List<pointIndexHit> hitInfo;
             labelList hitSurfaces;
             geometryToConformTo.findSurfaceNearest
             (
@@ -342,7 +630,6 @@ int main(int argc, char *argv[])
 
     mc.clean_all() ;
 
-
     Info<< "End\n" << endl;
 
     return 0;

applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/cvMeshSurfaceSimplify_non_octree.C

@@ -0,0 +1,352 @@
+/*---------------------------------------------------------------------------*\
+ =========                   |
+ \\      /   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
+    cvMeshSurfaceSimplify
+
+Description
+    Simplifies surfaces by resampling.
+
+    Uses Thomas Lewiner's topology preserving MarchingCubes.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "searchableSurfaces.H"
+#include "conformationSurfaces.H"
+#include "triSurfaceMesh.H"
+
+#include "MarchingCubes.h"
+
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Re-sample surfaces used in cvMesh operation"
+    );
+    //argList::validArgs.append("inputFile");
+    argList::validArgs.append("(nx ny nz)");
+    argList::validArgs.append("outputName");
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    runTime.functionObjects().off();
+
+    const Vector<label> n(IStringStream(args.args()[1])());
+    const fileName exportName = args.args()[2];
+
+    Info<< "Reading surfaces as specified in the cvMeshDict and"
+        << " writing re-sampled " << n << " to " << exportName
+        << nl << endl;
+
+    cpuTime timer;
+
+    IOdictionary cvMeshDict
+    (
+        IOobject
+        (
+            "cvMeshDict",
+            runTime.system(),
+            runTime,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    );
+
+    // Define/load all geometry
+    searchableSurfaces allGeometry
+    (
+        IOobject
+        (
+            "cvSearchableSurfaces",
+            runTime.constant(),
+            "triSurface",
+            runTime,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        ),
+        cvMeshDict.subDict("geometry")
+    );
+
+    Info<< "Geometry read in = "
+        << timer.cpuTimeIncrement() << " s." << nl << endl;
+
+
+    Random rndGen(64293*Pstream::myProcNo());
+
+    conformationSurfaces geometryToConformTo
+    (
+        runTime,
+        rndGen,
+        allGeometry,
+        cvMeshDict.subDict("surfaceConformation")
+    );
+
+    Info<< "Set up geometry in = "
+        << timer.cpuTimeIncrement() << " s." << nl << endl;
+
+
+
+    // Extend
+    treeBoundBox bb = geometryToConformTo.globalBounds();
+    {
+        const vector smallVec = 0.1*bb.span();
+        bb.min() -= smallVec;
+        bb.max() += smallVec;
+    }
+
+    Info<< "Meshing to bounding box " << bb << nl << endl;
+
+    const vector span(bb.span());
+    const vector d
+    (
+        span.x()/(n.x()-1),
+        span.y()/(n.y()-1),
+        span.z()/(n.z()-1)
+    );
+
+    MarchingCubes mc(span.x(), span.y(), span.z() ) ;
+    mc.set_resolution(n.x(), n.y(), n.z());
+    mc.init_all() ;
+
+
+    // Generate points
+    pointField points(mc.size_x()*mc.size_y()*mc.size_z());
+    label pointI = 0;
+
+    point pt;
+    for( int k = 0 ; k < mc.size_z() ; k++ )
+    {
+        pt.z() = bb.min().z() + k*d.z();
+        for( int j = 0 ; j < mc.size_y() ; j++ )
+        {
+            pt.y() = bb.min().y() + j*d.y();
+            for( int i = 0 ; i < mc.size_x() ; i++ )
+            {
+                pt.x() = bb.min().x() + i*d.x();
+                points[pointI++] = pt;
+            }
+        }
+    }
+
+    Info<< "Generated " << points.size() << " sampling points in = "
+        << timer.cpuTimeIncrement() << " s." << nl << endl;
+
+
+    // Compute data
+    const searchableSurfaces& geometry = geometryToConformTo.geometry();
+    const labelList& surfaces = geometryToConformTo.surfaces();
+
+    scalarField signedDist;
+    labelList nearestSurfaces;
+    searchableSurfacesQueries::signedDistance
+    (
+        geometry,
+        surfaces,
+        points,
+        scalarField(points.size(), sqr(GREAT)),
+        searchableSurface::OUTSIDE,     // for non-closed surfaces treat as
+                                        // outside
+        nearestSurfaces,
+        signedDist
+    );
+
+    // Fill elements
+    pointI = 0;
+    for( int k = 0 ; k < mc.size_z() ; k++ )
+    {
+        for( int j = 0 ; j < mc.size_y() ; j++ )
+        {
+            for( int i = 0 ; i < mc.size_x() ; i++ )
+            {
+                mc.set_data(float(signedDist[pointI++]), i, j, k);
+            }
+        }
+    }
+
+    Info<< "Determined signed distance in = "
+        << timer.cpuTimeIncrement() << " s." << nl << endl;
+
+
+    mc.run() ;
+
+    Info<< "Constructed iso surface in = "
+        << timer.cpuTimeIncrement() << " s." << nl << endl;
+
+
+    mc.clean_temps() ;
+
+
+
+    // Write output file
+    if (mc.ntrigs() > 0)
+    {
+        Triangle* triangles = mc.triangles();
+        List<labelledTri> tris(mc.ntrigs());
+        forAll(tris, triI)
+        {
+            tris[triI] = labelledTri
+            (
+                triangles[triI].v1,
+                triangles[triI].v2,
+                triangles[triI].v3,
+                0                       // region
+            );
+        }
+
+
+        Vertex* vertices = mc.vertices();
+        pointField points(mc.nverts());
+        forAll(points, pointI)
+        {
+            Vertex& v = vertices[pointI];
+            points[pointI] = point
+            (
+                bb.min().x() + v.x*span.x()/mc.size_x(),
+                bb.min().y() + v.y*span.y()/mc.size_y(),
+                bb.min().z() + v.z*span.z()/mc.size_z()
+            );
+        }
+
+
+        // Find correspondence to original surfaces
+        labelList regionOffsets(surfaces.size());
+        label nRegions = 0;
+        forAll(surfaces, i)
+        {
+            const wordList& regions = geometry[surfaces[i]].regions();
+            regionOffsets[i] = nRegions;
+            nRegions += regions.size();
+        }
+
+
+        geometricSurfacePatchList patches(nRegions);
+        nRegions = 0;
+        forAll(surfaces, i)
+        {
+            const wordList& regions = geometry[surfaces[i]].regions();
+
+            forAll(regions, regionI)
+            {
+                patches[nRegions] = geometricSurfacePatch
+                (
+                    "patch",
+                    geometry[surfaces[i]].name() + "_" + regions[regionI],
+                    nRegions
+                );
+                nRegions++;
+            }
+        }
+
+        triSurface s(tris, patches, points, true);
+
+        Info<< "Extracted triSurface in = "
+            << timer.cpuTimeIncrement() << " s." << nl << endl;
+
+
+        // Find out region on local surface of nearest point
+        {
+            List<pointIndexHit> hitInfo;
+            labelList hitSurfaces;
+            geometryToConformTo.findSurfaceNearest
+            (
+                s.faceCentres(),
+                scalarField(s.size(), sqr(GREAT)),
+                hitInfo,
+                hitSurfaces
+            );
+
+            // Get region
+            DynamicList<pointIndexHit> surfInfo(hitSurfaces.size());
+            DynamicList<label> surfIndices(hitSurfaces.size());
+
+            forAll(surfaces, surfI)
+            {
+                // Extract info on this surface
+                surfInfo.clear();
+                surfIndices.clear();
+                forAll(hitSurfaces, triI)
+                {
+                    if (hitSurfaces[triI] == surfI)
+                    {
+                        surfInfo.append(hitInfo[triI]);
+                        surfIndices.append(triI);
+                    }
+                }
+
+                // Calculate sideness of these surface points
+                labelList region;
+                geometry[surfaces[surfI]].getRegion(surfInfo, region);
+
+                forAll(region, i)
+                {
+                    label triI = surfIndices[i];
+                    s[triI].region() = regionOffsets[surfI]+region[i];
+                }
+            }
+        }
+
+        Info<< "Re-patched surface in = "
+            << timer.cpuTimeIncrement() << " s." << nl << endl;
+
+        triSurfaceMesh smesh
+        (
+            IOobject
+            (
+                exportName,
+                runTime.constant(), // instance
+                "triSurface",
+                runTime,            // registry
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            s
+        );
+
+        Info<< "writing surfMesh:\n  "
+            << smesh.searchableSurface::objectPath() << nl << endl;
+        smesh.searchableSurface::write();
+
+        Info<< "Written surface in = "
+            << timer.cpuTimeIncrement() << " s." << nl << endl;
+    }
+
+    mc.clean_all() ;
+
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -158,7 +158,8 @@ castellatedMeshControls
             //      internal : keep them as internal faces (default)
             //      baffle   : create baffles from them. This gives more
             //                 freedom in mesh motion
-            //      boundary : create loose-standing boundary faces.
+            //      boundary : create free-standing boundary faces (baffles
+            //                 but without the shared points)
             //faceType internal;
         }
     }

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C

@@ -525,12 +525,6 @@ int main(int argc, char *argv[])
             surfaceDict.lookupOrAddDefault<Switch>("writeVTK", "off");
         const Switch writeObj =
             surfaceDict.lookupOrAddDefault<Switch>("writeObj", "off");
-        const Switch writeFeatureEdgeMesh =
-            surfaceDict.lookupOrAddDefault<Switch>
-            (
-                "writeFeatureEdgeMesh",
-                "off"
-            );
 
         const Switch curvature =
             surfaceDict.lookupOrAddDefault<Switch>("curvature", "off");
@@ -773,6 +767,8 @@ int main(int argc, char *argv[])
         Info<< nl << "Writing extendedFeatureEdgeMesh to "
             << feMesh.objectPath() << endl;
 
+        mkDir(feMesh.path());
+
         if (writeObj)
         {
             feMesh.writeObj(feMesh.path()/surfFileName.lessExt().name());
@@ -781,29 +777,26 @@ int main(int argc, char *argv[])
         feMesh.write();
 
         // Write a featureEdgeMesh for backwards compatibility
-        if (writeFeatureEdgeMesh)
+        featureEdgeMesh bfeMesh
-        {
+        (
-            featureEdgeMesh bfeMesh
+            IOobject
             (
-                IOobject
+                surfFileName.lessExt().name() + ".eMesh",   // name
-                (
+                runTime.constant(),                         // instance
-                    surfFileName.lessExt().name() + ".eMesh",   // name
+                "triSurface",
-                    runTime.constant(),                         // instance
+                runTime,                                    // registry
-                    "triSurface",
+                IOobject::NO_READ,
-                    runTime,                                    // registry
+                IOobject::AUTO_WRITE,
-                    IOobject::NO_READ,
+                false
-                    IOobject::AUTO_WRITE,
+            ),
-                    false
+            feMesh.points(),
-                ),
+            feMesh.edges()
-                feMesh.points(),
+        );
-                feMesh.edges()
-            );
 
-            Info<< nl << "Writing featureEdgeMesh to "
+        Info<< nl << "Writing featureEdgeMesh to "
-                << bfeMesh.objectPath() << endl;
+            << bfeMesh.objectPath() << endl;
 
-            bfeMesh.regIOobject::write();
+        bfeMesh.regIOobject::write();
-        }
 
         triSurfaceMesh searchSurf
         (

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict

@@ -87,9 +87,6 @@ surface2.nas
 
     // Write options
 
-        // Write .eMesh file (for snappyHexMesh)
-        writeFeatureEdgeMesh    no;
-
         // Write features to obj format for postprocessing
         writeObj                yes;
 
@@ -99,5 +96,4 @@ surface2.nas
 }
 
 
-
 // ************************************************************************* //

applications/utilities/surface/surfaceMeshTriangulate/surfaceMeshTriangulate.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -312,7 +312,7 @@ int main(int argc, char *argv[])
 
             triSurface allSurf(allFaces, allPatches, allPoints);
 
-            // Cleanup (which does point merge as well
+            // Cleanup (which does point merge as well)
             allSurf.cleanup(false);
 
             // Write surface mesh

bin/tools/pre-commit-hook

@@ -275,6 +275,7 @@ checkCopyright()
     year=$(date +%Y)
     echo "$hookName: check copyright ..." 1>&2
 
+    badFiles=$(
     for f in $fileList
     do
         startYear=`grep "Copyright.*OpenFOAM" $f | sed 's/[^0-9]*\([0-9]*\).*/\1/g'`
@@ -287,19 +288,24 @@ checkCopyright()
                 # Date is of type 2011-2012 OpenFOAM Foundation
                 if [ "$year" != "$endYear" ]
                 then
-                    echo "Updated copyright for: $f"
+                    echo "Updated copyright for: $f" 1>&2
+                    echo "$f"
                     sed -i "s/$startYear-$endYear OpenFOAM/$startYear-$year OpenFOAM/g" $f
                 fi
             else
                 # Date is of type 2011 OpenFOAM Foundation
                 if [ "$year" != "$startYear" ]
                 then
-                    echo "Updated copyright for: $f"
+                    echo "$f"
+                    echo "Updated copyright for: $f" 1>&2
                     sed -i "s/$startYear OpenFOAM/$startYear-$year OpenFOAM/g" $f
                 fi
             fi
         fi
     done
+    )
+
+    dieOnBadFiles "Some copyright dates were automatically updated; Please check these before pushing"
 }
 
 

etc/config/settings.csh

@@ -41,6 +41,9 @@ alias _foamAddMan 'setenv MANPATH \!*\:${MANPATH}'
 # Set environment variables according to system type
 setenv WM_ARCH `uname -s`
 
+# Default WM_COMPILER_LIB_ARCH for 32bit
+setenv WM_COMPILER_LIB_ARCH
+
 switch ($WM_ARCH)
 case Linux:
     setenv WM_ARCH linux
@@ -269,12 +272,7 @@ case ThirdParty:
         _foamAddPath    $gccDir/bin
 
         # add compiler libraries to run-time environment
-        # 64-bit needs lib64, but 32-bit needs lib (not lib32)
+        _foamAddLib     $gccDir/lib$WM_COMPILER_LIB_ARCH
-        if ($WM_ARCH_OPTION == 64 && $?WM_COMPILER_LIB_ARCH) then
-            _foamAddLib     $gccDir/lib$WM_COMPILER_LIB_ARCH
-        else
-            _foamAddLib     $gccDir/lib
-        endif
 
         # add gmp/mpfr libraries to run-time environment
         _foamAddLib     $gmpDir/lib
@@ -371,21 +369,6 @@ unset boost_version cgal_version
 unsetenv MPI_ARCH_PATH MPI_HOME FOAM_MPI_LIBBIN
 
 switch ("$WM_MPLIB")
-case OPENMPI:
-    setenv FOAM_MPI openmpi-1.5.4
-    # optional configuration tweaks:
-    _foamSource `$WM_PROJECT_DIR/bin/foamEtcFile config/openmpi.csh`
-
-    setenv MPI_ARCH_PATH $WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$FOAM_MPI
-
-    # Tell OpenMPI where to find its install directory
-    setenv OPAL_PREFIX $MPI_ARCH_PATH
-
-    _foamAddPath    $MPI_ARCH_PATH/bin
-    _foamAddLib     $MPI_ARCH_PATH/lib
-    _foamAddMan     $MPI_ARCH_PATH/share/man
-    breaksw
-
 case SYSTEMOPENMPI:
     # Use the system installed openmpi, get library directory via mpicc
     setenv FOAM_MPI openmpi-system
@@ -410,13 +393,36 @@ case SYSTEMOPENMPI:
     unset libDir
     breaksw
 
+case OPENMPI:
+    setenv FOAM_MPI openmpi-1.5.4
+    # optional configuration tweaks:
+    _foamSource `$WM_PROJECT_DIR/bin/foamEtcFile config/openmpi.csh`
+
+    setenv MPI_ARCH_PATH $WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$FOAM_MPI
+
+    # Tell OpenMPI where to find its install directory
+    setenv OPAL_PREFIX $MPI_ARCH_PATH
+
+    _foamAddPath    $MPI_ARCH_PATH/bin
+
+    # 64-bit on OpenSuSE 12.1 uses lib64 others use lib
+    _foamAddLib     $MPI_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH
+    _foamAddLib     $MPI_ARCH_PATH/lib
+
+    _foamAddMan     $MPI_ARCH_PATH/share/man
+    breaksw
+
 case MPICH:
     setenv FOAM_MPI mpich2-1.1.1p1
     setenv MPI_HOME $WM_THIRD_PARTY_DIR/$FOAM_MPI
     setenv MPI_ARCH_PATH $WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$FOAM_MPI
 
     _foamAddPath    $MPI_ARCH_PATH/bin
+
+    # 64-bit on OpenSuSE 12.1 uses lib64 others use lib
+    _foamAddLib     $MPI_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH
     _foamAddLib     $MPI_ARCH_PATH/lib
+
     _foamAddMan     $MPI_ARCH_PATH/share/man
     breaksw
 
@@ -427,7 +433,11 @@ case MPICH-GM:
     setenv GM_LIB_PATH /opt/gm/lib64
 
     _foamAddPath    $MPI_ARCH_PATH/bin
+
+    # 64-bit on OpenSuSE 12.1 uses lib64 others use lib
+    _foamAddLib     $MPI_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH
     _foamAddLib     $MPI_ARCH_PATH/lib
+
     _foamAddLib     $GM_LIB_PATH
     breaksw
 

etc/config/settings.sh

@@ -294,14 +294,7 @@ OpenFOAM | ThirdParty)
         _foamAddPath    $gccDir/bin
 
         # add compiler libraries to run-time environment
-        # 64-bit needs lib64, but 32-bit needs lib (not lib32)
+        _foamAddLib     $gccDir/lib$WM_COMPILER_LIB_ARCH
-        if [ "$WM_ARCH_OPTION" = 64 ]
-        then
-            _foamAddLib     $gccDir/lib$WM_COMPILER_LIB_ARCH
-        else
-            _foamAddLib     $gccDir/lib
-        fi
-
 
         # add gmp/mpfr libraries to run-time environment
         _foamAddLib     $gmpDir/lib
@@ -401,21 +394,6 @@ unset boost_version cgal_version
 unset MPI_ARCH_PATH MPI_HOME FOAM_MPI_LIBBIN
 
 case "$WM_MPLIB" in
-OPENMPI)
-    export FOAM_MPI=openmpi-1.5.4
-    # optional configuration tweaks:
-    _foamSource `$WM_PROJECT_DIR/bin/foamEtcFile config/openmpi.sh`
-
-    export MPI_ARCH_PATH=$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$FOAM_MPI
-
-    # Tell OpenMPI where to find its install directory
-    export OPAL_PREFIX=$MPI_ARCH_PATH
-
-    _foamAddPath    $MPI_ARCH_PATH/bin
-    _foamAddLib     $MPI_ARCH_PATH/lib
-    _foamAddMan     $MPI_ARCH_PATH/share/man
-    ;;
-
 SYSTEMOPENMPI)
     # Use the system installed openmpi, get library directory via mpicc
     export FOAM_MPI=openmpi-system
@@ -441,13 +419,36 @@ SYSTEMOPENMPI)
     unset libDir
     ;;
 
+OPENMPI)
+    export FOAM_MPI=openmpi-1.5.4
+    # optional configuration tweaks:
+    _foamSource `$WM_PROJECT_DIR/bin/foamEtcFile config/openmpi.sh`
+
+    export MPI_ARCH_PATH=$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$FOAM_MPI
+
+    # Tell OpenMPI where to find its install directory
+    export OPAL_PREFIX=$MPI_ARCH_PATH
+
+    _foamAddPath    $MPI_ARCH_PATH/bin
+
+    # 64-bit on OpenSuSE 12.1 uses lib64 others use lib
+    _foamAddLib     $MPI_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH
+    _foamAddLib     $MPI_ARCH_PATH/lib
+
+    _foamAddMan     $MPI_ARCH_PATH/share/man
+    ;;
+
 MPICH)
     export FOAM_MPI=mpich2-1.1.1p1
     export MPI_HOME=$WM_THIRD_PARTY_DIR/$FOAM_MPI
     export MPI_ARCH_PATH=$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$FOAM_MPI
 
     _foamAddPath    $MPI_ARCH_PATH/bin
+
+    # 64-bit on OpenSuSE 12.1 uses lib64 others use lib
+    _foamAddLib     $MPI_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH
     _foamAddLib     $MPI_ARCH_PATH/lib
+
     _foamAddMan     $MPI_ARCH_PATH/share/man
     ;;
 
@@ -458,7 +459,11 @@ MPICH-GM)
     export GM_LIB_PATH=/opt/gm/lib64
 
     _foamAddPath    $MPI_ARCH_PATH/bin
+
+    # 64-bit on OpenSuSE 12.1 uses lib64 others use lib
+    _foamAddLib     $MPI_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH
     _foamAddLib     $MPI_ARCH_PATH/lib
+
     _foamAddLib     $GM_LIB_PATH
     ;;
 

src/OpenFOAM/db/IOstreams/StringStreams/OStringStream.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -116,14 +116,14 @@ public:
 
         // Edit
 
-            //- Clear the OStringStream
+            //- Rewind the OStringStream
             void rewind()
             {
-#               if __GNUC__ < 4 && __GNUC_MINOR__ < 4
+                #if __GNUC__ < 4 && __GNUC_MINOR__ < 4
                 stdStream().rdbuf()->pubsetbuf(" ", 1);
-#               else
+                #else
                 stdStream().rdbuf()->pubseekpos(0);
-#               endif
+                #endif
             }
 
 

src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.C

@@ -128,8 +128,6 @@ dimensioned<Type>::dimensioned
     dimensions_(dimSet),
     value_(pTraits<Type>::zero)
 {
-    Info<< "dimensioned<Type>::dimensioned" << endl;
-
     token nextToken(is);
     is.putBack(nextToken);
 

src/OpenFOAM/matrices/solution/solution.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,6 +65,8 @@ void Foam::solution::read(const dictionary& dict)
         else
         {
             // backwards compatibility
+            fieldRelaxDict_.clear();
+
             const wordList entryNames(relaxDict.toc());
             forAll(entryNames, i)
             {

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -544,7 +544,7 @@ public:
                 const labelListList& globalEdgeSlaves() const;
                 const labelListList& globalEdgeTransformedSlaves() const;
                 const mapDistribute& globalEdgeSlavesMap() const;
-                //- Is my edge same orientation master edge
+                //- Is my edge same orientation as master edge
                 const PackedBoolList& globalEdgeOrientation() const;
 
             // Coupled point to boundary faces. These are uncoupled boundary

src/OpenFOAM/meshes/polyMesh/polyMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1092,7 +1092,8 @@ const Foam::pointField& Foam::polyMesh::oldPoints() const
 
 Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
 (
-    const pointField& newPoints
+    const pointField& newPoints,
+    const bool deleteDemandDrivenData
 )
 {
     if (debug)
@@ -1146,6 +1147,14 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
 
     // Force recalculation of all geometric data with new points
 
+    if (deleteDemandDrivenData)
+    {
+        // Remove the stored tet base points
+        tetBasePtIsPtr_.clear();
+        // Remove the cell tree
+        cellTreePtr_.clear();
+    }
+
     bounds_ = boundBox(points_);
     boundary_.movePoints(points_);
 

src/OpenFOAM/meshes/polyMesh/polyMesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -442,7 +442,11 @@ public:
             }
 
             //- Move points, returns volumes swept by faces in motion
-            virtual tmp<scalarField> movePoints(const pointField&);
+            virtual tmp<scalarField> movePoints
+            (
+                const pointField&,
+                const bool deleteDemandDrivenData = false
+            );
 
             //- Reset motion
             void resetMotion() const;

src/OpenFOAM/meshes/polyMesh/polyMeshIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -407,9 +407,12 @@ Foam::polyMesh::readUpdateState Foam::polyMesh::readUpdate()
 
         clearGeom();
 
-        points_.instance() = pointsInst;
 
-        points_ = pointIOField
+        label nOldPoints = points_.size();
+
+        points_.clear();
+
+        pointIOField newPoints
         (
             IOobject
             (
@@ -423,6 +426,19 @@ Foam::polyMesh::readUpdateState Foam::polyMesh::readUpdate()
             )
         );
 
+        if (nOldPoints != 0 && nOldPoints != newPoints.size())
+        {
+            FatalErrorIn("polyMesh::readUpdate()")
+                << "Point motion detected but number of points "
+                << newPoints.size() << " in "
+                << newPoints.objectPath() << " does not correspond to "
+                << " current " << nOldPoints
+                << exit(FatalError);
+        }
+
+        points_.transfer(newPoints);
+        points_.instance() = pointsInst;
+
         // Derived info
         bounds_ = boundBox(points_);
 

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -197,7 +197,7 @@ Foam::label Foam::polyMeshTetDecomposition::findBasePoint
     }
 
     // If a base point hasn't triggered a return by now, then there is
-    // non that can produce a good decomposition
+    // none that can produce a good decomposition
     return -1;
 }
 

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,9 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    Calculation of shape function product for a tetrahedron
-
 \*---------------------------------------------------------------------------*/
 
 #include "tetrahedron.H"
@@ -32,6 +29,96 @@ Description
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Point, class PointRef>
+void Foam::tetrahedron<Point, PointRef>::tetOverlap
+(
+    const tetrahedron<Point, PointRef>& tetB,
+    tetIntersectionList& insideTets,
+    label& nInside,
+    tetIntersectionList& outsideTets,
+    label& nOutside
+) const
+{
+    // Work storage
+    tetIntersectionList cutInsideTets;
+    label nCutInside = 0;
+
+    nInside = 0;
+    storeOp inside(insideTets, nInside);
+    storeOp cutInside(cutInsideTets, nCutInside);
+
+    nOutside = 0;
+    storeOp outside(outsideTets, nOutside);
+
+
+    // Cut tetA with all inwards pointing faces of tetB. Any tets remaining
+    // in aboveTets are inside tetB.
+
+    {
+        // face0
+        plane pl0(tetB.b_, tetB.d_, tetB.c_);
+
+        // Cut and insert subtets into cutInsideTets (either by getting
+        // an index from freeSlots or by appending to insideTets) or
+        // insert into outsideTets
+        sliceWithPlane(pl0, cutInside, outside);
+    }
+
+    if (nCutInside == 0)
+    {
+        nInside = nCutInside;
+        return;
+    }
+
+    {
+        // face1
+        plane pl1(tetB.a_, tetB.c_, tetB.d_);
+
+        nInside = 0;
+
+        for (label i = 0; i < nCutInside; i++)
+        {
+            cutInsideTets[i].tet().sliceWithPlane(pl1, inside, outside);
+        }
+
+        if (nInside == 0)
+        {
+            return;
+        }
+    }
+
+    {
+        // face2
+        plane pl2(tetB.a_, tetB.d_, tetB.b_);
+
+        nCutInside = 0;
+
+        for (label i = 0; i < nInside; i++)
+        {
+            insideTets[i].tet().sliceWithPlane(pl2, cutInside, outside);
+        }
+
+        if (nCutInside == 0)
+        {
+            nInside = nCutInside;
+            return;
+        }
+    }
+
+    {
+        // face3
+        plane pl3(tetB.a_, tetB.b_, tetB.c_);
+
+        nInside = 0;
+
+        for (label i = 0; i < nCutInside; i++)
+        {
+            cutInsideTets[i].tet().sliceWithPlane(pl3, inside, outside);
+        }
+    }
+}
+
+
 // (Probably very inefficient) minimum containment sphere calculation.
 // From http://www.imr.sandia.gov/papers/imr11/shewchuk2.pdf:
 // Sphere ctr is smallest one of

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,6 +87,13 @@ class tetrahedron
 {
 public:
 
+    // Public typedefs
+
+        //- Storage type for tets originating from intersecting tets.
+        //  (can possibly be smaller than 200)
+        typedef FixedList<tetPoints, 200> tetIntersectionList;
+
+
         // Classes for use in sliceWithPlane. What to do with decomposition
         // of tet.
 
@@ -111,11 +118,11 @@ public:
             //- Store resulting tets
             class storeOp
             {
-                FixedList<tetPoints, 200>& tets_;
+                tetIntersectionList& tets_;
                 label& nTets_;
 
             public:
-                inline storeOp(FixedList<tetPoints, 200>&, label&);
+                inline storeOp(tetIntersectionList&, label&);
 
                 inline void operator()(const tetPoints&);
             };
@@ -261,6 +268,16 @@ public:
                 BelowTetOp& belowOp
             ) const;
 
+            //- Decompose tet into tets inside and outside other tet
+            inline void tetOverlap
+            (
+                const tetrahedron<Point, PointRef>& tetB,
+                tetIntersectionList& insideTets,
+                label& nInside,
+                tetIntersectionList& outsideTets,
+                label& nOutside
+            ) const;
+
 
             //- Return (min)containment sphere, i.e. the smallest sphere with
             //  all points inside. Returns pointHit with:

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H

@@ -556,7 +556,7 @@ inline void Foam::tetrahedron<Point, PointRef>::sumVolOp::operator()
 template<class Point, class PointRef>
 inline Foam::tetrahedron<Point, PointRef>::storeOp::storeOp
 (
-    FixedList<tetPoints, 200>& tets,
+    tetIntersectionList& tets,
     label& nTets
 )
 :

src/Pstream/mpi/UPstream.C

@@ -41,40 +41,6 @@ License
 #   define MPI_SCALAR MPI_DOUBLE
 #endif
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    typedef struct
-    {
-        scalar value;
-        label count;
-    } CountAndValue;
-
-    void reduceSum
-    (
-        void *in,
-        void *inOut,
-        int *len,
-        MPI_Datatype *dptr
-    )
-    {
-        CountAndValue* inPtr =
-            reinterpret_cast<CountAndValue*>(in);
-        CountAndValue* inOutPtr =
-            reinterpret_cast<CountAndValue*>(inOut);
-
-        for (int i=0; i< *len; ++i)
-        {
-            inOutPtr->value += inPtr->value;
-            inOutPtr->count += inPtr->count;
-            inPtr++;
-            inOutPtr++;
-        }
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // NOTE:
@@ -493,101 +459,11 @@ void Foam::sumReduce
     const int tag
 )
 {
-    static bool hasDataType_ = false;
+    vector2D twoScalars(Value, scalar(Count));
-    static MPI_Datatype mesg_mpi_strct_;
+    reduce(twoScalars, sumOp<vector2D>());
-    static MPI_Op myOp_;
 
-    if (Pstream::debug)
+    Value = twoScalars.x();
-    {
+    Count = twoScalars.y();
-        Pout<< "Foam::sumReduce : value:" << Value
-            << " count:" << Count << endl;
-    }
-
-    if (!UPstream::parRun())
-    {
-        return;
-    }
-
-    if (UPstream::nProcs() <= UPstream::nProcsSimpleSum)
-    {
-        reduce(Value, sumOp<scalar>(), tag);
-        reduce(Count, sumOp<label>(), tag);
-    }
-    else
-    {
-        CountAndValue in,out;
-
-        if (!hasDataType_)
-        {
-            int lengths[2];
-            lengths[0] = 1;
-            lengths[1] = 1;
-            MPI_Datatype types[2];
-            types[0] = MPI_DOUBLE;
-            types[1] = MPI_INT;
-            MPI_Aint addresses[2];
-            MPI_Address(&in.value, &addresses[0]);
-            MPI_Address(&in.count, &addresses[1]);
-            MPI_Aint offsets[2];
-            offsets[0] = 0;
-            offsets[1] = addresses[1]-addresses[0];
-
-            if
-            (
-                MPI_Type_create_struct
-                (
-                    2,
-                    lengths,
-                    offsets,
-                    types,
-                    &mesg_mpi_strct_
-                )
-            )
-            {
-                FatalErrorIn("Foam::sumReduce()")
-                    << "MPI_Type_create_struct" << abort(FatalError);
-            }
-            if (MPI_Type_commit(&mesg_mpi_strct_))
-            {
-                FatalErrorIn("Foam::sumReduce()")
-                    << "MPI_Type_commit" << abort(FatalError);
-            }
-            if (MPI_Op_create(reduceSum, true, &myOp_))
-            {
-                FatalErrorIn("Foam::sumReduce()")
-                    << "MPI_Op_create" << abort(FatalError);
-            }
-
-            hasDataType_ = true;
-        }
-
-        in.value = Value;
-        in.count = Count;
-        if
-        (
-            MPI_Allreduce
-            (
-                &in,
-                &out,
-                1,
-                mesg_mpi_strct_,
-                myOp_,
-                MPI_COMM_WORLD
-            )
-        )
-        {
-            FatalErrorIn("Foam::sumReduce(..)")
-                << "Problem." << abort(FatalError);
-        }
-        Value = out.value;
-        Count = out.count;
-    }
-
-    if (Pstream::debug)
-    {
-        Pout<< "Foam::reduce : reduced value:" << Value
-            << " reduced count:" << Count << endl;
-    }
 }
 
 

src/dynamicMesh/motionSmoother/motionSmoother.C

@@ -903,7 +903,9 @@ Foam::tmp<Foam::scalarField> Foam::motionSmoother::movePoints
         testSyncPositions(newPoints, 1E-6*mesh_.bounds().mag());
     }
 
-    tmp<scalarField> tsweptVol = mesh_.movePoints(newPoints);
+    // Move actual mesh points. Make sure to delete tetBasePtIs so it
+    // gets rebuilt.
+    tmp<scalarField> tsweptVol = mesh_.movePoints(newPoints, true);
 
     pp_.movePoints(mesh_.points());
 

src/fieldSources/Make/files

@@ -15,6 +15,10 @@ basicSource/rotorDiskSource/profileModel/profileModel.C
 basicSource/rotorDiskSource/profileModel/profileModelList.C
 basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C
 basicSource/rotorDiskSource/profileModel/series/seriesProfile.C
+basicSource/rotorDiskSource/trimModel/trimModel/trimModel.C
+basicSource/rotorDiskSource/trimModel/trimModel/trimModelNew.C
+basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C
+basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C
 
 basicSource/actuationDiskSource/actuationDiskSource.C
 basicSource/radialActuationDiskSource/radialActuationDiskSource.C
@@ -25,4 +29,4 @@ $(interRegion)/constantHeatTransfer/constantHeatTransfer.C
 $(interRegion)/tabulatedHeatTransfer/tabulatedHeatTransfer.C
 $(interRegion)/variableHeatTransfer/variableHeatTransfer.C
 
-LIB = $(FOAM_LIBBIN)/libfieldSources
+LIB = $(FOAM_LIBBIN)/libfieldSources

src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C

@@ -26,8 +26,8 @@ License
 #include "rotorDiskSource.H"
 #include "addToRunTimeSelectionTable.H"
 #include "mathematicalConstants.H"
+#include "trimModel.H"
 #include "unitConversion.H"
-#include "geometricOneField.H"
 #include "fvMatrices.H"
 #include "syncTools.H"
 
@@ -123,6 +123,8 @@ void Foam::rotorDiskSource::checkData()
 
 void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
 {
+    area_ = 0.0;
+
     static const scalar tol = 0.8;
 
     const label nInternalFaces = mesh_.nInternalFaces();
@@ -252,7 +254,7 @@ void Foam::rotorDiskSource::createCoordinateSystem()
     {
         case gmAuto:
         {
-            // determine rotation origin
+            // determine rotation origin (cell volume weighted)
             scalar sumV = 0.0;
             const scalarField& V = mesh_.V();
             const vectorField& C = mesh_.C();
@@ -262,6 +264,8 @@ void Foam::rotorDiskSource::createCoordinateSystem()
                 sumV += V[cellI];
                 origin += V[cellI]*C[cellI];
             }
+            reduce(origin, sumOp<vector>());
+            reduce(sumV, sumOp<scalar>());
             origin /= sumV;
 
             // determine first radial vector
@@ -277,6 +281,8 @@ void Foam::rotorDiskSource::createCoordinateSystem()
                     magR = mag(test);
                 }
             }
+            reduce(dx1, maxMagSqrOp<vector>());
+            magR = mag(dx1);
 
             // determine second radial vector and cross to determine axis
             forAll(cells_, i)
@@ -292,14 +298,18 @@ void Foam::rotorDiskSource::createCoordinateSystem()
                     }
                 }
             }
+            reduce(axis, maxMagSqrOp<vector>());
             axis /= mag(axis);
 
-            // axis direction is somewhat arbitrary - check if user needs
+            // correct the axis direction using a point above the rotor
-            // needs to reverse
-            bool reverse(readBool(coeffs_.lookup("reverseAxis")));
-            if (reverse)
             {
-                axis *= -1.0;
+                vector pointAbove(coeffs_.lookup("pointAbove"));
+                vector dir = pointAbove - origin;
+                dir /= mag(dir);
+                if ((dir & axis) < 0)
+                {
+                    axis *= -1.0;
+                }
             }
 
             coeffs_.lookup("refDirection") >> refDir;
@@ -375,9 +385,6 @@ void Foam::rotorDiskSource::constructGeometry()
         scalar cPos = cos(beta);
         scalar sPos = sin(beta);
         invR_[i] = tensor(cPos, 0.0, -sPos, 0.0, 1.0, 0.0, sPos, 0.0, cPos);
-
-        // geometric angle of attack - not including twist [radians]
-        alphag_[i] = trim_.alphaC + trim_.A*cos(psi) + trim_.B*sin(psi);
     }
 }
 
@@ -439,13 +446,12 @@ Foam::rotorDiskSource::rotorDiskSource
     inletVelocity_(vector::zero),
     tipEffect_(1.0),
     flap_(),
-    trim_(),
+    trim_(trimModel::New(*this, coeffs_)),
     blade_(coeffs_.subDict("blade")),
     profiles_(coeffs_.subDict("profiles")),
     x_(cells_.size(), vector::zero),
     R_(cells_.size(), I),
     invR_(cells_.size(), I),
-    alphag_(cells_.size(), 0.0),
     area_(cells_.size(), 0.0),
     coordSys_(false),
     rMax_(0.0)
@@ -462,34 +468,179 @@ Foam::rotorDiskSource::~rotorDiskSource()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::rotorDiskSource::calculate
+(
+    const vectorField& U,
+    const scalarField& alphag,
+    vectorField& force,
+    const bool divideVolume,
+    const bool output
+) const
+{
+    tmp<volScalarField> trho;
+    if (rhoName_ != "none")
+    {
+        trho = mesh_.lookupObject<volScalarField>(rhoName_);
+    }
+
+    const scalarField& V = mesh_.V();
+
+
+    // logging info
+    scalar dragEff = 0.0;
+    scalar liftEff = 0.0;
+    scalar AOAmin = GREAT;
+    scalar AOAmax = -GREAT;
+
+    forAll(cells_, i)
+    {
+        if (area_[i] > ROOTVSMALL)
+        {
+            const label cellI = cells_[i];
+
+            const scalar radius = x_[i].x();
+
+            // velocity in local cylindrical reference frame
+            vector Uc = coordSys_.localVector(U[cellI]);
+
+            // apply correction in local system due to coning
+            Uc = R_[i] & Uc;
+
+            // set radial component of velocity to zero
+            Uc.x() = 0.0;
+
+            // remove blade linear velocity from blade normal component
+            Uc.y() -= radius*omega_;
+
+            // determine blade data for this radius
+            // i1 = index of upper bound data point in blade list
+            scalar twist = 0.0;
+            scalar chord = 0.0;
+            label i1 = -1;
+            label i2 = -1;
+            scalar invDr = 0.0;
+            blade_.interpolate(radius, twist, chord, i1, i2, invDr);
+
+            // effective angle of attack
+            scalar alphaEff =
+                mathematical::pi + atan2(Uc.z(), Uc.y()) - (alphag[i] + twist);
+
+            if (alphaEff > mathematical::pi)
+            {
+                alphaEff -= mathematical::twoPi;
+            }
+            if (alphaEff < -mathematical::pi)
+            {
+                alphaEff += mathematical::twoPi;
+            }
+
+            AOAmin = min(AOAmin, alphaEff);
+            AOAmax = max(AOAmax, alphaEff);
+
+            // determine profile data for this radius and angle of attack
+            const label profile1 = blade_.profileID()[i1];
+            const label profile2 = blade_.profileID()[i2];
+
+            scalar Cd1 = 0.0;
+            scalar Cl1 = 0.0;
+            profiles_[profile1].Cdl(alphaEff, Cd1, Cl1);
+
+            scalar Cd2 = 0.0;
+            scalar Cl2 = 0.0;
+            profiles_[profile2].Cdl(alphaEff, Cd2, Cl2);
+
+            scalar Cd = invDr*(Cd2 - Cd1) + Cd1;
+            scalar Cl = invDr*(Cl2 - Cl1) + Cl1;
+
+            // apply tip effect for blade lift
+            scalar tipFactor = neg(radius/rMax_ - tipEffect_);
+
+            // calculate forces perpendicular to blade
+            scalar pDyn = 0.5*magSqr(Uc);
+            if (trho.valid())
+            {
+                pDyn *= trho()[cellI];
+            }
+
+            scalar f = pDyn*chord*nBlades_*area_[i]/mathematical::twoPi;
+            vector localForce = vector(0.0, f*Cd, tipFactor*f*Cl);
+
+            // convert force from local coning system into rotor cylindrical
+            localForce = invR_[i] & localForce;
+
+            // accumulate forces
+            dragEff += localForce.y();
+            liftEff += localForce.z();
+
+            // convert force to global cartesian co-ordinate system
+            force[cellI] = coordSys_.globalVector(localForce);
+
+            if (divideVolume)
+            {
+                force[cellI] /= V[cellI];
+            }
+        }
+    }
+
+
+    if (output)
+    {
+        reduce(AOAmin, minOp<scalar>());
+        reduce(AOAmax, maxOp<scalar>());
+        reduce(dragEff, sumOp<scalar>());
+        reduce(liftEff, sumOp<scalar>());
+
+        Info<< type() << " output:" << nl
+            << "    min/max(AOA)   = " << radToDeg(AOAmin) << ", "
+            << radToDeg(AOAmax) << nl
+            << "    Effective drag = " << dragEff << nl
+            << "    Effective lift = " << liftEff << endl;
+    }
+}
+
+
 void Foam::rotorDiskSource::addSup(fvMatrix<vector>& eqn, const label fieldI)
 {
-    // add source to rhs of eqn
+    dimensionSet dims = dimless;
-
-    const volVectorField& U = eqn.psi();
-
     if (eqn.dimensions() == dimForce)
     {
         coeffs_.lookup("rhoName") >> rhoName_;
-
+        dims.reset(dimForce/dimVolume);
-        const volScalarField& rho =
-            mesh_.lookupObject<volScalarField>(rhoName_);
-
-        eqn -= calculateForces
-            (
-                rho.internalField(),
-                inflowVelocity(U),
-                dimForce/dimVolume
-            );
     }
     else
     {
-        eqn -= calculateForces
+        dims.reset(dimForce/dimVolume/dimDensity);
-            (
+    }
-                oneField(),
+
-                inflowVelocity(U),
+    volVectorField force
-                dimForce/dimVolume/dimDensity
+    (
-            );
+        IOobject
+        (
+            "rotorForce",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector("zero", dims, vector::zero)
+    );
+
+    const volVectorField& U = eqn.psi();
+
+    const vectorField Uin = inflowVelocity(U);
+
+    trim_->correct(Uin, force);
+
+    calculate(Uin, trim_->thetag(), force);
+
+
+    // add source to rhs of eqn
+    eqn -= force;
+
+    if (mesh_.time().outputTime())
+    {
+        force.write();
     }
 }
 
@@ -522,12 +673,10 @@ bool Foam::rotorDiskSource::read(const dictionary& dict)
         flapCoeffs.lookup("beta1") >> flap_.beta1;
         flapCoeffs.lookup("beta2") >> flap_.beta2;
         flap_.beta0 = degToRad(flap_.beta0);
+        flap_.beta1 = degToRad(flap_.beta1);
+        flap_.beta2 = degToRad(flap_.beta2);
 
-        const dictionary& trimCoeffs(coeffs_.subDict("trimCoeffs"));
+        trim_->read(coeffs_);        
-        trimCoeffs.lookup("alphaC") >> trim_.alphaC;
-        trimCoeffs.lookup("A") >> trim_.A;
-        trimCoeffs.lookup("B") >> trim_.B;
-        trim_.alphaC = degToRad(trim_.alphaC);
 
         checkData();
 
@@ -537,7 +686,7 @@ bool Foam::rotorDiskSource::read(const dictionary& dict)
 
         if (debug)
         {
-            writeField("alphag", alphag_, true);
+            writeField("alphag", trim_->thetag()(), true);
             writeField("faceArea", area_, true);
         }
 

src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,13 +37,16 @@ Description
             fieldNames      (U);    // names of fields on which to apply source
             rhoName         rho;    // density field if compressible case
             nBlades         3;      // number of blades
-            tip effect      0.96;   // normalised radius above which lift = 0
+            tipEffect       0.96;   // normalised radius above which lift = 0
 
             inletFlowType   local;  // inlet flow type specification
 
             geometryMode    auto;   // geometry specification
 
             refDirection    (-1 0 0); // reference direction
+                                      // - used as reference for psi angle
+
+            trimModel       fixed;  // fixed || targetForce
 
             flapCoeffs
             {
@@ -51,16 +54,12 @@ Description
                 beta1           0;  // lateral flapping coeff
                 beta2           0;  // longitudinal flapping coeff
             }
-            trimCoeffs
+
-            {
-                alphac          15; // collective pitch angle [deg]
-                A               0;  // lateral cyclic coeff
-                B               0;  // longitudinal cyclic coeff
-            }
             blade
             {
                 ...
             }
+
             profiles
             {
                 ...
@@ -102,6 +101,9 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
+class trimModel;
+
 /*---------------------------------------------------------------------------*\
                       Class rotorDiskSource Declaration
 \*---------------------------------------------------------------------------*/
@@ -140,13 +142,6 @@ protected:
             scalar beta2;   // longitudinal flapping coeff
         };
 
-        struct trimData
-        {
-            scalar alphaC;  // collective pitch angle
-            scalar A;       // lateral cyclic coeff
-            scalar B;       // longitudinal cyclic coeff
-        };
-
 
     // Protected data
 
@@ -172,8 +167,8 @@ protected:
         //- Blade flap coefficients [rad/s]
         flapData flap_;
 
-        //- Blad trim coefficients
+        //- Trim model
-        trimData trim_;
+        autoPtr<trimModel> trim_;
 
         //- Blade data
         bladeModel blade_;
@@ -181,7 +176,8 @@ protected:
         //- Profile data
         profileModelList profiles_;
 
-        //- Cell centre positions in local rotor frame (Cartesian x, y, z)
+        //- Cell centre positions in local rotor frame
+        //  (Cylindrical r, theta, z)
         List<point> x_;
 
         //- Rotation tensor for flap angle
@@ -190,9 +186,6 @@ protected:
         //- Inverse rotation tensor for flap angle
         List<tensor> invR_;
 
-        //- Geometric angle of attack [deg]
-        List<scalar> alphag_;
-
         //- Area [m2]
         List<scalar> area_;
 
@@ -220,15 +213,6 @@ protected:
         //- Return the inlet flow field
         tmp<vectorField> inflowVelocity(const volVectorField& U) const;
 
-        //- Calculate forces
-        template<class RhoType>
-        tmp<volVectorField> calculateForces
-        (
-            const RhoType& rho,
-            const vectorField& U,
-            const dimensionSet& dims
-        );
-
         //- Helper function to write rotor values
         template<class Type>
         void writeField
@@ -264,6 +248,29 @@ public:
 
     // Member Functions
 
+        // Access
+            
+            //- Return the cell centre positions in local rotor frame
+            // (Cylindrical r, theta, z)
+            inline const List<point>& x() const;
+
+            //- Return the rotor co-ordinate system (r, theta, z)
+            inline const cylindricalCS& coordSys() const;
+
+
+        // Evaluation
+
+            //- Calculate forces
+            void calculate
+            (
+                const vectorField& U,
+                const scalarField& alphag,
+                vectorField& force,
+                const bool divideVolume = true,
+                const bool output = true
+            ) const;
+
+
         // Source term addition
 
             //- Source term to fvMatrix<vector>
@@ -286,6 +293,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "rotorDiskSourceI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #ifdef NoRepository
     #include "rotorDiskSourceTemplates.C"
 #endif

src/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceI.H

@@ -0,0 +1,41 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "rotorDiskSource.H"
+
+const Foam::List<Foam::point>& Foam::rotorDiskSource::x() const
+{
+    return x_;
+}
+
+
+const Foam::cylindricalCS& Foam::rotorDiskSource::coordSys() const
+{
+    return coordSys_;
+}
+
+
+// ************************************************************************* //
+

src/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceTemplates.C

@@ -24,13 +24,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "rotorDiskSource.H"
-#include "addToRunTimeSelectionTable.H"
-#include "mathematicalConstants.H"
-#include "unitConversion.H"
 #include "volFields.H"
 
-using namespace Foam::constant::mathematical;
-
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class Type>
@@ -81,149 +76,4 @@ void Foam::rotorDiskSource::writeField
 }
 
 
-template<class RhoType>
-Foam::tmp<Foam::volVectorField> Foam::rotorDiskSource::calculateForces
-(
-    const RhoType& rho,
-    const vectorField& U,
-    const dimensionSet& dims
-)
-{
-    tmp<volVectorField> tForce
-    (
-        new volVectorField
-        (
-            IOobject
-            (
-                "rotorForce",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedVector("zero", dims, vector::zero)
-        )
-    );
-
-    vectorField& force = tForce().internalField();
-    const scalarField& V = mesh_.V();
-
-
-    // logging info
-    scalar dragEff = 0.0;
-    scalar liftEff = 0.0;
-    scalar AOAmin = GREAT;
-    scalar AOAmax = -GREAT;
-
-    forAll(cells_, i)
-    {
-        if (area_[i] > ROOTVSMALL)
-        {
-            const label cellI = cells_[i];
-
-            const scalar radius = x_[i].x();
-
-            // velocity in local cylindrical reference frame
-            vector Uc = coordSys_.localVector(U[cellI]);
-
-            // apply correction in local system due to coning
-            Uc = R_[i] & Uc;
-
-            // set radial component of velocity to zero
-            Uc.x() = 0.0;
-
-            // remove blade linear velocity from blade normal component
-            Uc.y() -= radius*omega_;
-
-            // velocity magnitude
-            scalar magUc = mag(Uc);
-
-            // determine blade data for this radius
-            // i1 = index of upper bound data point in blade list
-            scalar twist = 0.0;
-            scalar chord = 0.0;
-            label i1 = -1;
-            label i2 = -1;
-            scalar invDr = 0.0;
-            blade_.interpolate(radius, twist, chord, i1, i2, invDr);
-
-            // effective angle of attack
-            scalar alphaEff = pi + atan2(Uc.z(), Uc.y()) - (alphag_[i] + twist);
-            if (alphaEff > pi)
-            {
-                alphaEff -= twoPi;
-            }
-            if (alphaEff < -pi)
-            {
-                alphaEff += twoPi;
-            }
-
-            AOAmin = min(AOAmin, alphaEff);
-            AOAmax = max(AOAmax, alphaEff);
-
-            // determine profile data for this radius and angle of attack
-            const label profile1 = blade_.profileID()[i1];
-            const label profile2 = blade_.profileID()[i2];
-
-            scalar Cd1 = 0.0;
-            scalar Cl1 = 0.0;
-            profiles_[profile1].Cdl(alphaEff, Cd1, Cl1);
-
-            scalar Cd2 = 0.0;
-            scalar Cl2 = 0.0;
-            profiles_[profile2].Cdl(alphaEff, Cd2, Cl2);
-
-            scalar Cd = invDr*(Cd2 - Cd1) + Cd1;
-            scalar Cl = invDr*(Cl2 - Cl1) + Cl1;
-
-            // apply tip effect for blade lift
-            scalar tipFactor = 1.0;
-            if (radius/rMax_ > tipEffect_)
-            {
-                tipFactor = 0.0;
-            }
-
-            // calculate forces perpendicular to blade
-            scalar pDyn = 0.5*rho[cellI]*sqr(magUc);
-            scalar f = pDyn*chord*nBlades_*area_[i]/twoPi;
-            vector localForce = vector(0.0, f*Cd, tipFactor*f*Cl);
-
-            // convert force from local coning system into rotor cylindrical
-            localForce = invR_[i] & localForce;
-
-            // accumulate forces
-            dragEff += localForce.y();
-            liftEff += localForce.z();
-
-            // convert force to global cartesian co-ordinate system
-            force[cellI] = coordSys_.globalVector(localForce);
-
-            force[cellI] /= V[cellI];
-        }
-    }
-
-
-    if (mesh_.time().outputTime())
-    {
-        tForce().write();
-    }
-
-
-    reduce(AOAmin, minOp<scalar>());
-    reduce(AOAmax, maxOp<scalar>());
-    reduce(dragEff, sumOp<scalar>());
-    reduce(liftEff, sumOp<scalar>());
-
-    Info<< type() << " output:" << nl
-        << "    min/max(AOA)   = " << radToDeg(AOAmin) << ", "
-        << radToDeg(AOAmax) << nl
-        << "    Effective drag = " << dragEff << nl
-        << "    Effective lift = " << liftEff << endl;
-
-
-    return tForce;
-}
-
-
 // ************************************************************************* //

src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C

@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fixedTrim.H"
+#include "addToRunTimeSelectionTable.H"
+#include "unitConversion.H"
+#include "mathematicalConstants.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(fixedTrim, 0);
+
+    addToRunTimeSelectionTable(trimModel, fixedTrim, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fixedTrim::fixedTrim(const rotorDiskSource& rotor, const dictionary& dict)
+:
+    trimModel(rotor, dict, typeName),
+    thetag_(rotor.cells().size(), 0.0)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
+
+Foam::fixedTrim::~fixedTrim()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fixedTrim::read(const dictionary& dict)
+{
+    trimModel::read(dict);
+
+    scalar theta0 = degToRad(readScalar(coeffs_.lookup("theta0")));
+    scalar theta1c = degToRad(readScalar(coeffs_.lookup("theta1c")));
+    scalar theta1s = degToRad(readScalar(coeffs_.lookup("theta1s")));
+
+    const List<vector>& x = rotor_.x();
+    forAll(thetag_, i)
+    {
+        scalar psi = x[i].y();
+        if (psi < 0)
+        {
+            psi += mathematical::twoPi;
+        }
+
+        thetag_[i] = theta0 + theta1c*cos(psi) + theta1s*sin(psi);
+    }
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::fixedTrim::thetag() const
+{
+    return tmp<scalarField>(thetag_);
+}
+
+
+void Foam::fixedTrim::correct(const vectorField& U, vectorField& force)
+{
+    // do nothing
+}
+
+
+// ************************************************************************* //

src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.H

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::fixedTrim
+
+Description
+    Fixed trim coefficients
+
+SourceFiles
+    fixedTrim.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedTrim_H
+#define fixedTrim_H
+
+#include "trimModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class fixedTrim Declaration
+\*---------------------------------------------------------------------------*/
+
+class fixedTrim
+:
+    public trimModel
+{
+
+protected:
+
+    // Protected data
+
+        //- Geometric angle of attack [rad]
+        scalarField thetag_;
+
+
+public:
+
+    //- Run-time type information
+    TypeName("fixedTrim");
+
+    //- Constructor
+    fixedTrim(const rotorDiskSource& rotor, const dictionary& dict);
+
+    //- Destructor
+    virtual ~fixedTrim();
+
+
+    // Member functions
+
+        //- Read
+        void read(const dictionary& dict);
+
+        //- Return the geometric angle of attack [rad]
+        virtual tmp<scalarField> thetag() const;
+
+        //- Correct the model
+        virtual void correct(const vectorField& U, vectorField& force);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C

@@ -0,0 +1,234 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "targetForceTrim.H"
+#include "addToRunTimeSelectionTable.H"
+#include "unitConversion.H"
+#include "mathematicalConstants.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(targetForceTrim, 0);
+
+    addToRunTimeSelectionTable(trimModel, targetForceTrim, dictionary);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::vector Foam::targetForceTrim::calcForce
+(
+    const vectorField& U,
+    const scalarField& thetag,
+    vectorField& force    
+) const
+{
+    rotor_.calculate(U, thetag, force, false, false);
+
+    const labelList& cells = rotor_.cells();
+    const vectorField& C = rotor_.mesh().C();
+
+    const vector& origin = rotor_.coordSys().origin();
+    const vector& rollAxis = rotor_.coordSys().e1();
+    const vector& pitchAxis = rotor_.coordSys().e2();
+    const vector& yawAxis = rotor_.coordSys().e3();
+
+    vector f(vector::zero);
+    forAll(cells, i)
+    {
+        label cellI = cells[i];
+
+        vector moment = force[cellI]^(C[cellI] - origin);
+        f[0] += force[cellI] & yawAxis;
+        f[1] += moment & pitchAxis;
+        f[2] += moment & rollAxis;
+    }
+
+    reduce(f, sumOp<vector>());
+    
+    return f;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::targetForceTrim::targetForceTrim
+(
+    const rotorDiskSource& rotor,
+    const dictionary& dict
+)
+:
+    trimModel(rotor, dict, typeName),
+    calcFrequency_(-1),
+    target_(vector::zero),
+    theta_(vector::zero),
+    nIter_(50),
+    tol_(1e-8),
+    relax_(1.0),
+    dTheta_(degToRad(0.1))
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
+
+Foam::targetForceTrim::~targetForceTrim()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::targetForceTrim::read(const dictionary& dict)
+{
+    trimModel::read(dict);
+
+    const dictionary& targetDict(coeffs_.subDict("target"));
+    target_[0] = readScalar(targetDict.lookup("fThrust"));
+    target_[1] = readScalar(targetDict.lookup("mPitch"));
+    target_[2] = readScalar(targetDict.lookup("mRoll"));
+
+    const dictionary& pitchAngleDict(coeffs_.subDict("pitchAngles"));
+    theta_[0] = degToRad(readScalar(pitchAngleDict.lookup("theta0Ini")));
+    theta_[1] = degToRad(readScalar(pitchAngleDict.lookup("theta1cIni")));
+    theta_[2] = degToRad(readScalar(pitchAngleDict.lookup("theta1sIni")));
+
+    coeffs_.lookup("calcFrequency") >> calcFrequency_;
+
+    coeffs_.readIfPresent("nIter", nIter_);
+    coeffs_.readIfPresent("tol", tol_);
+    coeffs_.readIfPresent("relax", relax_);
+
+    if (coeffs_.readIfPresent("dTheta", dTheta_))
+    {
+        dTheta_ = degToRad(dTheta_);
+    }
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::targetForceTrim::thetag() const
+{
+    const List<vector>& x = rotor_.x();
+
+    tmp<scalarField> ttheta(new scalarField(x.size()));
+    scalarField& t = ttheta();
+
+    forAll(t, i)
+    {
+        scalar psi = x[i].y();
+        if (psi < 0)
+        {
+            psi += mathematical::twoPi;
+        }
+
+        t[i] = theta_[0] + theta_[1]*cos(psi) + theta_[2]*sin(psi);
+    }
+
+    return ttheta;
+}
+
+
+void Foam::targetForceTrim::correct(const vectorField& U, vectorField& force)
+{
+    if (rotor_.mesh().time().timeIndex() % calcFrequency_ == 0)
+    {
+        // iterate to find new pitch angles to achieve target force
+        scalar err = GREAT;
+        label iter = 0;
+        tensor J(tensor::zero);
+
+        while ((err > tol_) && (iter < nIter_))
+        {
+            // cache initial theta vector
+            vector theta0(theta_);
+
+            // set initial values
+            vector old = calcForce(U, thetag(), force);
+
+            // construct Jacobian by perturbing the pitch angles
+            // by +/-(dTheta_/2)
+            for (label pitchI = 0; pitchI < 3; pitchI++)
+            {
+                theta_[pitchI] -= dTheta_/2.0;
+                vector f0 = calcForce(U, thetag(), force);
+
+                theta_[pitchI] += dTheta_;
+                vector f1 = calcForce(U, thetag(), force);
+
+                vector ddTheta = (f1 - f0)/dTheta_;
+
+                J[pitchI + 0] = ddTheta[0];
+                J[pitchI + 3] = ddTheta[1];
+                J[pitchI + 6] = ddTheta[2];
+
+                theta_ = theta0;
+            }
+
+            // calculate the change in pitch angle vector
+            vector dt = inv(J) & (target_ - old);
+
+            // update pitch angles
+            vector thetaNew = theta_ + relax_*dt;
+
+            // update error
+            err = mag(thetaNew - theta_);
+
+            // update for next iteration
+            theta_ = thetaNew;
+            iter++;
+        }
+
+        if (iter == nIter_)
+        {
+            WarningIn
+            (
+                "void Foam::targetForceTrim::correct"
+                "("
+                    "const vectorField&, "
+                    "vectorField&"
+                ")"
+            )   << "Trim routine not converged in " << iter
+                << " iterations, max residual = " << err << endl;
+        }
+        else
+        {
+            Info<< type() << ": converged in " << iter
+                << " iterations" << endl;
+        }
+
+        Info<< "    new pitch angles:" << nl
+            << "        theta0  = " << radToDeg(theta_[0]) << nl
+            << "        theta1c = " << radToDeg(theta_[1]) << nl
+            << "        theta1s = " << radToDeg(theta_[2]) << nl
+            << endl;
+    }
+}
+
+
+// ************************************************************************* //

src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.H

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::targetForceTrim
+
+Description
+    Target force trim coefficients
+
+SourceFiles
+    targetForceTrim.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef targetForceTrim_H
+#define targetForceTrim_H
+
+#include "trimModel.H"
+#include "tensor.H"
+#include "vector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class targetForceTrim Declaration
+\*---------------------------------------------------------------------------*/
+
+class targetForceTrim
+:
+    public trimModel
+{
+
+protected:
+
+    // Protected data
+
+        //- Number of iterations between calls to 'correct'
+        label calcFrequency_;
+
+        //- Target force [N]
+        vector target_;
+
+        //- Pitch angles (collective, roll, pitch) [rad]
+        vector theta_;
+
+        //- Maximum number of iterations in trim routine
+        label nIter_;
+
+        //- Convergence tolerance
+        scalar tol_;
+
+        //- Under-relaxation coefficient
+        scalar relax_;
+
+        //- Perturbation angle used to determine jacobian
+        scalar dTheta_;
+
+
+    // Protected member functions
+
+        //- Calculate the rotor forces
+        vector calcForce
+        (
+            const vectorField& U,
+            const scalarField& alphag,
+            vectorField& force    
+        ) const;
+
+
+public:
+
+    //- Run-time type information
+    TypeName("targetForceTrim");
+
+    //- Constructor
+    targetForceTrim(const rotorDiskSource& rotor, const dictionary& dict);
+
+    //- Destructor
+    virtual ~targetForceTrim();
+
+
+    // Member functions
+
+        //- Read
+        void read(const dictionary& dict);
+
+        //- Return the geometric angle of attack [rad]
+        virtual tmp<scalarField> thetag() const;
+
+        //- Correct the model
+        virtual void correct(const vectorField& U, vectorField& force);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fieldSources/basicSource/rotorDiskSource/trimModel/trimModel/trimModel.C

@@ -0,0 +1,67 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "trimModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(trimModel, 0);
+    defineRunTimeSelectionTable(trimModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::trimModel::trimModel
+(
+    const rotorDiskSource& rotor,
+    const dictionary& dict,
+    const word& name
+)
+:
+    rotor_(rotor),
+    name_(name),
+    coeffs_(dictionary::null)
+{
+    read(dict);
+}
+
+// * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
+
+Foam::trimModel::~trimModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::trimModel::read(const dictionary& dict)
+{
+    coeffs_ = dict.subDict(name_ + "Coeffs");
+}
+
+
+// ************************************************************************* //

src/fieldSources/basicSource/rotorDiskSource/trimModel/trimModel/trimModel.H

@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::trimModel
+
+Description
+    Trim model base class
+
+SourceFiles
+    trimModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef trimModel_H
+#define trimModel_H
+
+#include "rotorDiskSource.H"
+#include "dictionary.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class trimModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class trimModel
+{
+
+protected:
+
+    // Protected data
+
+        //- Reference to the rotor source model
+        const rotorDiskSource& rotor_;
+
+        //- Name of model
+        const word name_;
+
+        //- Coefficients dictionary
+        dictionary coeffs_;
+
+
+public:
+
+    //- Run-time type information
+    TypeName("trimModel");
+
+
+    // Declare runtime constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            trimModel,
+            dictionary,
+            (
+                const rotorDiskSource& rotor,
+                const dictionary& dict
+            ),
+            (rotor, dict)
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        trimModel
+        (
+            const rotorDiskSource& rotor,
+            const dictionary& dict,
+            const word& name
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected trim model
+        static autoPtr<trimModel> New
+        (
+            const rotorDiskSource& rotor,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~trimModel();
+
+
+    // Member functions
+
+        //- Read
+        virtual void read(const dictionary& dict);
+
+        //- Return the geometric angle of attack [rad]
+        virtual tmp<scalarField> thetag() const = 0;
+
+        //- Correct the model
+        virtual void correct(const vectorField& U, vectorField& force) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fieldSources/basicSource/rotorDiskSource/trimModel/trimModel/trimModelNew.C

@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "trimModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::trimModel> Foam::trimModel::New
+(
+    const rotorDiskSource& rotor,
+    const dictionary& dict
+)
+{
+    const word modelType(dict.lookup(typeName));
+
+    Info<< "    Selecting " << typeName << " " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "trimModel::New(const rotorDiskSource&, const dictionary&)"
+        )   << "Unknown " << typeName << " type "
+            << modelType << nl << nl
+            << "Valid " << typeName << " types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<trimModel>(cstrIter()(rotor, dict));
+}
+
+
+// ************************************************************************* //

src/fileFormats/vtk/vtkUnstructuredReader.C

@@ -36,14 +36,18 @@ defineTypeNameAndDebug(Foam::vtkUnstructuredReader, 0);
 
 template<>
 const char*
-Foam::NamedEnum<Foam::vtkUnstructuredReader::vtkDataType, 4>::names[] =
+Foam::NamedEnum<Foam::vtkUnstructuredReader::vtkDataType, 8>::names[] =
 {
     "int",
+    "unsigned_int",
+    "long",
+    "unsigned_long",
     "float",
+    "double",
     "string",
     "vtkIdType"
 };
-const Foam::NamedEnum<Foam::vtkUnstructuredReader::vtkDataType, 4>
+const Foam::NamedEnum<Foam::vtkUnstructuredReader::vtkDataType, 8>
 Foam::vtkUnstructuredReader::vtkDataTypeNames;
 
 
@@ -385,6 +389,9 @@ void Foam::vtkUnstructuredReader::readField
     switch (vtkDataTypeNames[dataType])
     {
         case VTK_INT:
+        case VTK_UINT:
+        case VTK_LONG:
+        case VTK_ULONG:
         case VTK_ID:
         {
             autoPtr<labelIOField> fieldVals
@@ -406,6 +413,7 @@ void Foam::vtkUnstructuredReader::readField
         break;
 
         case VTK_FLOAT:
+        case VTK_DOUBLE:
         {
             autoPtr<scalarIOField> fieldVals
             (
@@ -627,7 +635,7 @@ void Foam::vtkUnstructuredReader::read(ISstream& inFile)
             }
 
             word primitiveTag(inFile);
-            if (primitiveTag != "float")
+            if (primitiveTag != "float" && primitiveTag != "double")
             {
                 FatalIOErrorIn("vtkUnstructuredReader::read(..)", inFile)
                     << "Expected 'float' entry but found "
@@ -809,7 +817,11 @@ void Foam::vtkUnstructuredReader::read(ISstream& inFile)
                 3*wantedSize
             );
 
-            if (vtkDataTypeNames[dataType] == VTK_FLOAT)
+            if
+            (
+                vtkDataTypeNames[dataType] == VTK_FLOAT
+             || vtkDataTypeNames[dataType] == VTK_DOUBLE
+            )
             {
                 objectRegistry::iterator iter = reg.find(dataName);
                 scalarField s(*dynamic_cast<const scalarField*>(iter()));

src/fileFormats/vtk/vtkUnstructuredReader.H

@@ -28,6 +28,8 @@ Description
     Reader for vtk unstructured_grid legacy files. Supports single CELLS, POINTS
     etc. entry only.
 
+    - all integer types (int, unsigned_int, long etc.) become Foam::label
+    - all real types (float, double) become Foam::scalar
     - POINTS becomes OpenFOAM points
     - CELLS gets split into OpenFOAM
         - cells
@@ -69,12 +71,16 @@ public:
         enum vtkDataType
         {
             VTK_INT,
+            VTK_UINT,
+            VTK_LONG,
+            VTK_ULONG,
             VTK_FLOAT,
+            VTK_DOUBLE,
             VTK_STRING,
             VTK_ID
         };
 
-        static const NamedEnum<vtkDataType, 4> vtkDataTypeNames;
+        static const NamedEnum<vtkDataType, 8> vtkDataTypeNames;
 
 
         //- Enumeration defining the vtk dataset types

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -2216,28 +2216,6 @@ void Foam::meshRefinement::dumpRefinementLevel() const
 
         pointRefLevel.write();
     }
-
-    // Dump cell centres
-    {
-        for (direction i=0; i<vector::nComponents; i++)
-        {
-            volScalarField cci
-            (
-                IOobject
-                (
-                    "cc" + word(vector::componentNames[i]),
-                    mesh_.time().timeName(),
-                    mesh_,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    false
-                ),
-                mesh_.C().component(i)
-            );
-
-            cci.write();
-        }
-    }
 }
 
 

src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshWeights.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -118,7 +118,7 @@ void Foam::meshToMesh::calculateInverseVolumeWeights() const
     inverseVolumeWeightsPtr_ = new scalarListList(toMesh_.nCells());
     scalarListList& invVolCoeffs = *inverseVolumeWeightsPtr_;
 
-    labelListList& cellToCell = *cellToCellAddressingPtr_;
+    const labelListList& cellToCell = cellToCellAddressing();
 
     tetOverlapVolume overlapEngine;
 
@@ -129,8 +129,8 @@ void Foam::meshToMesh::calculateInverseVolumeWeights() const
         if (overlapCells.size() > 0)
         {
             invVolCoeffs[celli].setSize(overlapCells.size());
-            scalar v(0);
+
-            forAll (overlapCells, j)
+            forAll(overlapCells, j)
             {
                 label cellFrom = overlapCells[j];
                 treeBoundBox bbFromMesh
@@ -142,7 +142,7 @@ void Foam::meshToMesh::calculateInverseVolumeWeights() const
                     )
                 );
 
-                v = overlapEngine.cellCellOverlapVolumeMinDecomp
+                scalar v = overlapEngine.cellCellOverlapVolumeMinDecomp
                 (
                     toMesh_,
                     celli,
@@ -151,19 +151,14 @@ void Foam::meshToMesh::calculateInverseVolumeWeights() const
                     cellFrom,
                     bbFromMesh
                 );
-                invVolCoeffs[celli][j] =  v/toMesh_.V()[celli];
+                invVolCoeffs[celli][j] = v/toMesh_.V()[celli];
-            }
-            if (celli == 2)
-            {
-                Info << "cellToCell :" << cellToCell[celli] << endl;
-                Info << "invVolCoeffs :" << invVolCoeffs[celli] << endl;
             }
         }
     }
 }
 
 
-void  Foam::meshToMesh::calculateCellToCellAddressing() const
+void Foam::meshToMesh::calculateCellToCellAddressing() const
 {
     if (debug)
     {

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C

@@ -45,102 +45,15 @@ Foam::tetOverlapVolume::tetOverlapVolume()
 
 // * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
 
-void Foam::tetOverlapVolume::tetTetOverlap
-(
-    const tetPoints& tetA,
-    const tetPoints& tetB,
-    FixedList<tetPoints, 200>& insideTets,
-    label& nInside,
-    FixedList<tetPoints, 200>& outsideTets,
-    label& nOutside
-) const
-{
-    // Work storage
-    FixedList<tetPoints, 200> cutInsideTets;
-    label nCutInside = 0;
-
-    tetPointRef::storeOp inside(insideTets, nInside);
-    tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
-    tetPointRef::dummyOp outside;
-
-
-    // Cut tetA with all inwards pointing faces of tetB. Any tets remaining
-    // in aboveTets are inside tetB.
-
-    {
-        // face0
-        plane pl0(tetB[1], tetB[3], tetB[2]);
-
-        // Cut and insert subtets into cutInsideTets (either by getting
-        // an index from freeSlots or by appending to insideTets) or
-        // insert into outsideTets
-        tetA.tet().sliceWithPlane(pl0, cutInside, outside);
-    }
-
-    if (nCutInside == 0)
-    {
-        nInside = nCutInside;
-        return;
-    }
-
-    {
-        // face1
-        plane pl1(tetB[0], tetB[2], tetB[3]);
-
-        nInside = 0;
-
-        for (label i = 0; i < nCutInside; i++)
-        {
-            cutInsideTets[i].tet().sliceWithPlane(pl1, inside, outside);
-        }
-
-        if (nInside == 0)
-        {
-            return;
-        }
-    }
-
-    {
-        // face2
-        plane pl2(tetB[0], tetB[3], tetB[1]);
-
-        nCutInside = 0;
-
-        for (label i = 0; i < nInside; i++)
-        {
-            insideTets[i].tet().sliceWithPlane(pl2, cutInside, outside);
-        }
-
-        if (nCutInside == 0)
-        {
-            nInside = nCutInside;
-            return;
-        }
-    }
-
-    {
-        // face3
-        plane pl3(tetB[0], tetB[1], tetB[2]);
-
-        nInside = 0;
-
-        for (label i = 0; i < nCutInside; i++)
-        {
-            cutInsideTets[i].tet().sliceWithPlane(pl3, inside, outside);
-        }
-    }
-}
-
-
 Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol
 (
     const tetPoints& tetA,
     const tetPoints& tetB
 ) const
 {
-    FixedList<tetPoints, 200> insideTets;
+    tetPointRef::tetIntersectionList insideTets;
     label nInside = 0;
-    FixedList<tetPoints, 200> cutInsideTets;
+    tetPointRef::tetIntersectionList cutInsideTets;
     label nCutInside = 0;
 
     tetPointRef::storeOp inside(insideTets, nInside);
@@ -222,7 +135,7 @@ Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp
     const primitiveMesh& meshB,
     const label cellBI,
     const treeBoundBox& cellBbB
-)
+) const
 {
     const cell& cFacesA = meshA.cells()[cellAI];
     const point& ccA = meshA.cellCentres()[cellAI];

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H

@@ -55,17 +55,6 @@ class tetOverlapVolume
 {
     // Private member functions
 
-        //- Tet overlap
-        void tetTetOverlap
-        (
-            const tetPoints& tetA,
-            const tetPoints& tetB,
-            FixedList<tetPoints, 200>& insideTets,
-            label& nInside,
-            FixedList<tetPoints, 200>& outsideTets,
-            label& nOutside
-        ) const;
-
         //- Tet Overlap Vol
         scalar tetTetOverlapVol
         (
@@ -115,7 +104,7 @@ public:
             const primitiveMesh& meshB,
             const label cellBI,
             const treeBoundBox& cellBbB
-        );
+        ) const;
 };
 
 

src/turbulenceModels/compressible/RAS/LRR/LRR.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -351,6 +351,7 @@ void LRR::correct()
     (
         fvm::ddt(rho_, epsilon_)
       + fvm::div(phi_, epsilon_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
     //- fvm::laplacian(Ceps*rho_*(k_/epsilon_)*R_, epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
@@ -393,6 +394,7 @@ void LRR::correct()
     (
         fvm::ddt(rho_, R_)
       + fvm::div(phi_, R_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), R_)
     //- fvm::laplacian(Cs*rho_*(k_/epsilon_)*R_, R_)
       - fvm::laplacian(DREff(), R_)
       + fvm::Sp(Clrr1_*rho_*epsilon_/k_, R_)

src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -389,6 +389,7 @@ void LaunderGibsonRSTM::correct()
     (
         fvm::ddt(rho_, epsilon_)
       + fvm::div(phi_, epsilon_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
     //- fvm::laplacian(Ceps*rho_*(k_/epsilon_)*R_, epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
@@ -432,6 +433,7 @@ void LaunderGibsonRSTM::correct()
     (
         fvm::ddt(rho_, R_)
       + fvm::div(phi_, R_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), R_)
     //- fvm::laplacian(Cs*rho_*(k_/epsilon_)*R_, R_)
       - fvm::laplacian(DREff(), R_)
       + fvm::Sp(Clg1_*rho_*epsilon_/k_, R_)

src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -314,6 +314,7 @@ void LaunderSharmaKE::correct()
     (
         fvm::ddt(rho_, epsilon_)
       + fvm::div(phi_, epsilon_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
         C1_*G*epsilon_/k_ + fvm::SuSp((C3_ - 2.0/3.0*C1_)*rho_*divU, epsilon_)
@@ -333,6 +334,7 @@ void LaunderSharmaKE::correct()
     (
         fvm::ddt(rho_, k_)
       + fvm::div(phi_, k_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(), k_)
      ==
         G - fvm::SuSp(2.0/3.0*rho_*divU, k_)

src/turbulenceModels/compressible/RAS/RNGkEpsilon/RNGkEpsilon.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -321,6 +321,7 @@ void RNGkEpsilon::correct()
     (
         fvm::ddt(rho_, epsilon_)
       + fvm::div(phi_, epsilon_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
       ==
         (C1_ - R)*G*epsilon_/k_
@@ -342,6 +343,7 @@ void RNGkEpsilon::correct()
     (
         fvm::ddt(rho_, k_)
       + fvm::div(phi_, k_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(), k_)
       ==
         G - fvm::SuSp(2.0/3.0*rho_*divU, k_)

src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.C

@@ -416,6 +416,7 @@ void SpalartAllmaras::correct()
     (
         fvm::ddt(rho_, nuTilda_)
       + fvm::div(phi_, nuTilda_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), nuTilda_)
       - fvm::laplacian(DnuTildaEff(), nuTilda_)
       - Cb2_/sigmaNut_*rho_*magSqr(fvc::grad(nuTilda_))
      ==

src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -292,6 +292,7 @@ void kEpsilon::correct()
     (
         fvm::ddt(rho_, epsilon_)
       + fvm::div(phi_, epsilon_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
         C1_*G*epsilon_/k_
@@ -313,6 +314,7 @@ void kEpsilon::correct()
     (
         fvm::ddt(rho_, k_)
       + fvm::div(phi_, k_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(), k_)
      ==
         G

src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -411,6 +411,7 @@ void kOmegaSST::correct()
     (
         fvm::ddt(rho_, omega_)
       + fvm::div(phi_, omega_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), omega_)
       - fvm::laplacian(DomegaEff(F1), omega_)
      ==
         rhoGammaF1*GbyMu
@@ -435,6 +436,7 @@ void kOmegaSST::correct()
     (
         fvm::ddt(rho_, k_)
       + fvm::div(phi_, k_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(F1), k_)
      ==
         min(G, (c1_*betaStar_)*rho_*k_*omega_)

src/turbulenceModels/compressible/RAS/realizableKE/realizableKE.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -331,6 +331,7 @@ void realizableKE::correct()
     (
         fvm::ddt(rho_, epsilon_)
       + fvm::div(phi_, epsilon_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
         C1*rho_*magS*epsilon_
@@ -355,6 +356,7 @@ void realizableKE::correct()
     (
         fvm::ddt(rho_, k_)
       + fvm::div(phi_, k_)
+      - fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), k_)
       - fvm::laplacian(DkEff(), k_)
      ==
         G - fvm::SuSp(2.0/3.0*rho_*divU, k_)

tutorials/compressible/rhoSimplecFoam/squareBend/0/p

@@ -26,8 +26,8 @@ boundaryField
     }
     inlet
     {
-        //type            zeroGradient;
+        type            zeroGradient;
-        type            mixed;
+        //type            mixed;
         refValue        uniform 110000;
         refGradient     uniform 0;
         valueFraction   uniform 0.3;

tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSolution

@@ -65,8 +65,9 @@ relaxationFactors
     }
     equations
     {
+        p               1;
         U               0.9;
-        h               0.8;
+        h               0.9;
         k               0.9;
         epsilon         0.9;
     }

tutorials/multiphase/compressibleTwoPhaseEulerFoam/mixerVessel2D/system/topoSetDict

@@ -0,0 +1,32 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      topoSetDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    {
+        name    rotor;
+        type    cellSet;
+        action  new;
+        source  zoneToCell;
+        sourceInfo
+        {
+            name rotor;
+        }
+    }
+);
+
+// ************************************************************************* //

tutorials/multiphase/twoPhaseEulerFoam/mixerVessel2D/system/topoSetDict

@@ -0,0 +1,32 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      topoSetDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    {
+        name    rotor;
+        type    cellSet;
+        action  new;
+        source  zoneToCell;
+        sourceInfo
+        {
+            name rotor;
+        }
+    }
+);
+
+// ************************************************************************* //

wmake/rules/General/mplibMPICH

@@ -1,3 +1,3 @@
 PFLAGS     = -DMPICH_SKIP_MPICXX
 PINC       = -I$(MPI_ARCH_PATH)/include
-PLIBS      = -L$(MPI_ARCH_PATH)/lib -lmpich -lrt
+PLIBS      = -L$(MPI_ARCH_PATH)/lib$(WM_COMPILER_LIB_ARCH) -L$(MPI_ARCH_PATH)/lib -lmpich -lrt

wmake/rules/General/mplibMPICH-GM

@@ -1,3 +1,3 @@
 PFLAGS     =
 PINC       = -I$(MPI_ARCH_PATH)/include
-PLIBS      = -L$(MPI_ARCH_PATH)/lib -lmpich -L$(GM_LIB_PATH) -lgm
+PLIBS      = -L$(MPI_ARCH_PATH)/lib$(WM_COMPILER_LIB_ARCH) -L$(MPI_ARCH_PATH)/lib -lmpich -L$(GM_LIB_PATH) -lgm

wmake/rules/General/mplibOPENMPI

@@ -1,3 +1,3 @@
 PFLAGS     = -DOMPI_SKIP_MPICXX
 PINC       = -I$(MPI_ARCH_PATH)/include
-PLIBS      = -L$(MPI_ARCH_PATH)/lib -lmpi
+PLIBS      = -L$(MPI_ARCH_PATH)/lib$(WM_COMPILER_LIB_ARCH) -L$(MPI_ARCH_PATH)/lib -lmpi