[Rojj]

I hope someone can provide some advice...


I am meshing a large city geometry with snappyHexMesh and it keeps getting stuck at

Splitting mesh at surface intersections
---------------------------------------

Introducing baffles for 18354934 faces that are intersected by the surface.

There is no error and the CPUs are all running at 100%.



Honestly, I am a bit lost and I am not sure what to check.

I have been waiting for around an hour. The coarser version of the mesh for the same geometry goes through this part in a few minutes.

If helpful, this is the snappyHexMeshDict

// -*- C++ -*-
// File generated by PyFoam - sorry for the ugliness

FoamFile
{
 version 2.0;
 format ascii;
 class dictionary;
 object snappyHexMeshDict;
}

castellatedMesh yes;
snap yes;
addLayers no;
geometry
{
  refinement_1H
  {
    type searchableCylinder;
    point1 (-10.712 413.668 -109.993);
    point2 (-10.712 413.668 42.916);
    radius 204.29628787216436;
  }
  refinement_3H
  {
    type searchableCylinder;
    point1 (-10.712 413.668 -109.993);
    point2 (-10.712 413.668 37.9068);
    radius 306.4444318082466;
  }
  refinement_5H
  {
    type searchableCylinder;
    point1 (-10.712 413.668 -109.993);
    point2 (-10.712 413.668 32.9068);
    radius 408.5925757443287;
  }
  refinement_regions
  {
    type triSurfaceMesh;
    file "refinement_regions.stl";
  }
  mct-terrain.stl
  {
    name mct-terrain;
    type triSurfaceMesh;
  }
  fixed_surroundings.stl
  {
    name fixed_surroundings;
    type triSurfaceMesh;
  }
  far_buildings.stl
  {
    name far_buildings;
    type triSurfaceMesh;
  }
  boundary.stl
  {
    name boundary;
    regions
    {
      Bottom
      {
        name Bottom;
      }
      InletOutlet_0
      {
        name InletOutlet_0;
      }
      InletOutlet_1
      {
        name InletOutlet_1;
      }
      InletOutlet_2
      {
        name InletOutlet_2;
      }
      InletOutlet_3
      {
        name InletOutlet_3;
      }
      InletOutlet_4
      {
        name InletOutlet_4;
      }
      InletOutlet_5
      {
        name InletOutlet_5;
      }
      InletOutlet_6
      {
        name InletOutlet_6;
      }
      InletOutlet_7
      {
        name InletOutlet_7;
      }
      Top
      {
        name Top;
      }
    }
    type triSurfaceMesh;
  }
}     // Settings for the castellatedMesh generation.

castellatedMeshControls
{
  maxLocalCells 3859376;
  maxGlobalCells 247000000;     // Refinement parameters// ~~~~~~~~~~~~~~~~~~~~~// If local number of cells is >= maxLocalCells on any processor// switches from from refinement followed by balancing// (current method) to (weighted) balancing before refinement.// Overall cell limit (approximately). Refinement will stop immediately// upon reaching this number so a refinement level might not complete.// Note that this is the number of cells before removing the part which// is not 'visible' from the keepPoint. The final number of cells might// actually be a lot less.  // The surface refinement loop might spend lots of iterations refining just a// few cells. This setting will cause refinement to stop if <= minimumRefine// are selected for refinement. Note: it will at least do one iteration// (unless the number of cells to refine is 0)
  minRefinementCells 4;     // Allow a certain level of imbalance during refining// (since balancing is quite expensive)// Expressed as fraction of perfect balance (= overall number of cells /// nProcs). 0=balance always.
  maxLoadUnbalance 0.1;     // Number of buffer layers between different levels.// 1 means normal 2:1 refinement restriction, larger means slower// refinement.
  nCellsBetweenLevels 1;     // Explicit feature edge refinement// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~// Specifies a level for any cell intersected by explicitly provided// edges.// This is a featureEdgeMesh, read from constant/triSurface for now.// Specify 'levels' in the same way as the 'distance' mode in the// refinementRegions (see below). The old specification//      level   2;// is equivalent to//      levels  ((0 2));
  features
    (
    );     // Surface based refinement// ~~~~~~~~~~~~~~~~~~~~~~~~// Specifies two levels for every surface. The first is the minimum level,// every cell intersecting a surface gets refined up to the minimum level.// The second level is the maximum level. Cells that 'see' multiple// intersections where the intersections make an// angle > resolveFeatureAngle get refined up to the maximum level.// Surface-wise min and max refinement level
  refinementSurfaces
  {
    mct-terrain
    {
      level
        (
          2
          4
        );
    }
    boundary
    {
      level
        (
          2
          4
        );
    }
    fixed_surroundings
    {
      level
        (
          2
          3
        );
    }
    far_buildings
    {
      level
        (
          1
          2
        );
    }     //- Optional increment (on top of max level) in small gaps//gapLevelIncrement 2;//- Optional angle to detect small-large cell situation//  perpendicular to the surface. Is the angle of face w.r.t.//  the local surface normal. Use on flat(ish) surfaces only.//  Otherwise leave out or set to negative number.//perpendicularAngle 10;//- Optional faceZone and (for closed surface) cellZone with//  how to select the cells that are in the cellZone//  (inside / outside / specified insidePoint)//  The orientation of the faceZone is//  - if on cellZone(s) : point out of (maximum) cellZone//  - if freestanding   : oriented according to surface//faceZone sphere;//cellZone sphere;//cellZoneInside inside;  //outside/insidePoint//- Optional specification of what to do with faceZone faces://      internal : keep them as internal faces (default)//      baffle   : create baffles from them. This gives more//                 freedom in mesh motion//      boundary : create free-standing boundary faces (baffles//                 but without the shared points)//faceType baffle;
  }     // Feature angle:// - used if min and max refinement level of a surface differ// - used if feature snapping (see snapControls below) is used
  resolveFeatureAngle 30;     //- Optional increment (on top of max level) in small gaps//gapLevelIncrement 2;// Planar angle:// - used to determine if surface normals//   are roughly the same or opposite. Used//      - in proximity refinement//      - to decide when to merge free-standing baffles//        (if e.g. running in surfaceSimplify mode set this to 180 to//         merge all baffles)//      - in snapping to avoid snapping to nearest on 'wrong' side//        of thin gap//// If not specified same as resolveFeatureAngle
  planarAngle 30;     // Region-wise refinement// ~~~~~~~~~~~~~~~~~~~~~~// Specifies refinement level for cells in relation to a surface. One of// three modes// - distance. 'levels' specifies per distance to the surface the//   wanted refinement level. The distances need to be specified in//   increasing order.// - inside. 'levels' is only one entry and only the level is used. All//   cells inside the surface get refined up to the level. The surface//   needs to be closed for this to be possible.// - outside. Same but cells outside.
  refinementRegions
  {
    refinement_1H
    {
      mode inside;
      levels
        (

          (
            1000000000000000.0
            3
          )
        );
    }
    refinement_3H
    {
      mode inside;
      levels
        (

          (
            1000000000000000.0
            2
          )
        );
    }
    refinement_5H
    {
      mode inside;
      levels
        (

          (
            1000000000000000.0
            1
          )
        );
    }
    mct-terrain
    {
      mode distance;
      levels
        (

          (
            2.0
            3
          )
        );
    }
    fixed_surroundings
    {
      mode distance;
      levels
        (

          (
            1.0
            3
          )

          (
            4.0
            2
          )
        );
    }
    far_buildings
    {
      mode distance;
      levels
        (

          (
            1.0
            1
          )

          (
            4.0
            1
          )
        );
    }
  }     // Mesh selection// ~~~~~~~~~~~~~~// After refinement patches get added for all refinementSurfaces and// all cells intersecting the surfaces get put into these patches. The// section reachable from the locationInMesh is kept.// NOTE: This point should never be on a face, always inside a cell, even// after refinement.
  locationInMesh (-10.712045025920519 413.66830664447565 929.9970092773438);     // Whether any faceZones (as specified in the refinementSurfaces)// are only on the boundary of corresponding cellZones or also allow// free-standing zone faces. Not used if there are no faceZones.
  allowFreeStandingZoneFaces no;
}     // Settings for the snapping.

snapControls
{
  nSmoothPatch 3;
  tolerance 2.0;     // Number of patch smoothing iterations before finding correspondence// to surface// Maximum relative distance for points to be attracted by surface.// True distance is this factor times local maximum edge length.// Note: changed(corrected) w.r.t 17x! (17x used 2* tolerance)  // Number of mesh displacement relaxation iterations.
  nSolveIter 300;     // Maximum number of snapping relaxation iterations. Should stop// before upon reaching a correct mesh.
  nRelaxIter 10;     // Feature snapping// Number of feature edge snapping iterations.// Leave out altogether to disable.
  nFeatureSnapIter 15;     // Detect (geometric only) features by sampling the surface// (default=false).
  implicitFeatureSnap yes;
  explicitFeatureSnap no;
  multiRegionFeatureSnap no;
}     // Settings for the layer addition.

addLayersControls
{
  relativeSizes no;
  expansionRatio 1.05;     // Are the thickness parameters below relative to the undistorted// size of the refined cell outside layer (true) or absolute sizes (false).// Layer thickness specification. This can be specified in one of following// ways:// - expansionRatio and finalLayerThickness (cell nearest internal mesh)// - expansionRatio and firstLayerThickness (cell on surface)// - overall thickness and firstLayerThickness// - overall thickness and finalLayerThickness// - overall thickness and expansionRatio//// Note: the mode thus selected is global, i.e. one cannot override the//       mode on a per-patch basis (only the values can be overridden)// Expansion factor for layer mesh// Wanted thickness of the layer furthest away from the wall.// If relativeSizes this is relative to undistorted size of cell// outside layer.
  finalLayerThickness 0.25;     // Wanted overall thickness of layers.// If relativeSizes this is relative to undistorted size of cell// outside layer.//thickness 0.5// Minimum overall thickness of total layers. If for any reason layer// cannot be above minThickness do not add layer.// If relativeSizes this is relative to undistorted size of cell// outside layer..
  minThickness 0.1;     // Per final patch (so not geometry!) the layer information// Note: This behaviour changed after 21x. Any non-mentioned patches//       now slide unless://          - nSurfaceLayers is explicitly mentioned to be 0.//          - angle to nearest surface < slipFeatureAngle (see below)
  layers
  {
    mct-terrain
    {
      nSurfaceLayers 4;
    }
    fixed_surroundings
    {
      nSurfaceLayers 0;
    }
    far_buildings
    {
      nSurfaceLayers 0;
    }
  }     // If points get not extruded do nGrow layers of connected faces that are// also not grown. This helps convergence of the layer addition process// close to features.// Note: changed(corrected) w.r.t 17x! (didn't do anything in 17x)
  nGrow 0;     // Advanced settings// Static analysis of starting mesh// When not to extrude surface. 0 is flat surface, 90 is when two faces// are perpendicular
  featureAngle 310;     // Stop layer growth on highly warped cells
  maxFaceThicknessRatio 1.0;     // Patch displacement// Number of smoothing iterations of surface normals
  nSmoothSurfaceNormals 1;     // Smooth layer thickness over surface patches
  nSmoothThickness 10;     // Medial axis analysis// Angle used to pick up medial axis points// Note: changed(corrected) w.r.t 17x! 90 degrees corresponds to 130// in 17x.
  minMedianAxisAngle 5;     // Reduce layer growth where ratio thickness to medial// distance is large
  maxThicknessToMedialRatio 0.3;     // Number of smoothing iterations of interior mesh movement direction
  nSmoothNormals 3;     // Optional: limit the number of steps walking away from the surface.// Default is unlimited.//nMedialAxisIter 10;// Optional: smooth displacement after medial axis determination.// default is 0.//nSmoothDisplacement 90;// (wip)Optional: do not extrude a point if none of the surrounding points is// not extruded. Default is false.//detectExtrusionIsland true;        // Mesh shrinking        // Optional: at non-patched sides allow mesh to slip if extrusion// direction makes angle larger than slipFeatureAngle. Default is// 0.5*featureAngle.
  slipFeatureAngle 80;     // Maximum number of snapping relaxation iterations. Should stop// before upon reaching a correct mesh.
  nRelaxIter 3;     // Create buffer region for new layer terminations
  nBufferCellsNoExtrude 0;     // Overall max number of layer addition iterations. The mesher will// exit if it reaches this number of iterations; possibly with an// illegal mesh.
  nLayerIter 10;     // Max number of iterations after which relaxed meshQuality controls// get used. Up to nRelaxedIter it uses the settings in// meshQualityControls,// after nRelaxedIter it uses the values in// meshQualityControls::relaxed.
  nRelaxedIter 20;     // Additional reporting: if there are just a few faces where there// are mesh errors (after adding the layers) print their face centres.// This helps in tracking down problematic mesh areas.//additionalReporting true;
}     // Generic mesh quality settings. At any undoable phase these determine// where to undo.

meshQualityControls
{
#include "meshQualityDict"
  relaxed
  {
    maxNonOrtho 75;
  }     // Specify mesh quality constraints in separate dictionary so can// be reused (e.g. checkMesh -meshQuality)// Optional : some meshing phases allow usage of relaxed rules.// See e.g. addLayersControls::nRelaxedIter.// Maximum non-orthogonality allowed. Set to 180 to disable.// Advanced// Number of error distribution iterations
  nSmoothScale 4;     // amount to scale back displacement at error points
  errorReduction 0.75;
}     // Advanced//// Debug flags//debugFlags//(//    mesh            // write intermediate meshes//    intersections   // write current mesh intersections as .obj files//    featureSeeds    // write information about explicit feature edge//                    // refinement//    attraction      // write attraction as .obj files//    layerInfo       // write information about layers//);//

writeFlags
  (
    scalarLevels
    layerSets
    layerFields
  );     // write volScalarField with cellLevel for postprocessing// write cellSets, faceSets of faces in layer// write volScalarField for layer coverage// Merge tolerance. Is fraction of overall bounding box of initial mesh.// Note: the write tolerance needs to be higher than this.

mergeTolerance 1e-06;     // ************************************************************************* //

EDIT
The "problem" seems to go away by setting the levels in the fixed_surrounding refinement surfaces to a value of (3 3) rather than (2 3).

    fixed_surroundings
    {
      level
        (
          3 (instead of 2)
          3
        );
    }

I made this change based on a suggestion of a colleague, but I do not really understand why it works.

[Rojj]

Just in case it helps. It seems that improving the quality of the STL triSurfaces also solves the problem.

I am using MeshLab to re-triangulate. The filter I use is called Remeshing: Isotropic Explicit Remeshing