Incorrect geometry creation – forward-facing 2D step

[vthk]

Dear Developers,

I wish to report a bug that I noticed in the creation of the step geometry for 2D flows. I have tried this with both serial and mpi versions but the problem remains. The description of the problem is as follows:

When I try to create a step geometry for a 200-by-20 lattice units (lu) simulation domain with the 10 lu high step specified at 50 lu along the x-axis, I find the following:

1. If I pick the numOfCuboids = 1, I just get a channel and no step.
2. If I pick the numOfCuboids = 4, I get the step exactly at 50 lu.
3. If I pick the numOfCuboids = 3, I get the step at 66.67 lu.
3. If I pick the numOfCuboids = 5, I get the step at exactly 80 lu.

The same happens for the mpi-version, if one just replaces the number of cuboids by the number of threads/procs used with mpirun. All simulations run properly and complete.

I have tried this on a windows-10 laptop in the powershell (Ubuntu-20.04) and also on a linux workstation with Ubuntu LTS 18.04. For mpi, I am using openMpi and the compiler is g++ or mpic++.

If you need any more information for reproducing the problem at your end please do let me know. I would also very much appreciate it if you can let me know of a work-around for this.

The problem appears to resolve itself for the 3D-step geometry but I still need to test that more exhaustively.

Thanks.

Best regards,
User-vthk.

To allow for an easy reproduction of the problem, my code is posted below:
************************************************************************************
************************************************************************************
/* fstep2d.cpp:
* The implementation of a backward facing step. It is furthermore
* shown how to use checkpointing to save the state of the
* simulation regularly.
*/

#include “olb2D.h”
#ifndef OLB_PRECOMPILED // Unless precompiled version is used,
#include “olb2D.hh” // include full template code
#endif

#include <cmath>
#include <vector>
#include <iostream>
#include <fstream>
#include <iomanip>

using namespace olb;
using namespace olb::descriptors;

typedef double T;
//#define MRT
//#ifdef MRT
//#define DESCRIPTOR ForcedMRTD2Q9Descriptor
//#else
//#define DESCRIPTOR D2Q9<FORCE>
//#endif

#define DESCRIPTOR D2Q9<>

//typedef enum {forced, nonForced} FlowType;

// Partial Slip does not work with the superGeometry based creation of simulation domain,
// probably need to create boundary conditions and geometry explicitly. Even local BC is also
// unstable, use with Caution!
typedef enum {bounceBack, local, interpolated, freeSlip, partialSlip} BoundaryType; //
typedef enum {constantVel, constantPress, timeVaryingVel} InletType;

BoundaryType boundaryType = interpolated; //
InletType inletType = constantPress; // constantVel
//FlowType flowType = forced;

const T lx1 = 50.0e-9; // length of step in meter
const T ly1 = 10.0e-9; // height of step in meter
const T lx0 = 200.0e-9; // length of channel in meter
const T ly0 = 20.0e-9; // height of channel in meter
const int N = 10; // resolution of the model

const T dx = lx1/N;
const T dt = 5e-12;
const T density = 808.5;
const T uMax = 0.2;
const T outletPressure = 0; // in N/m^2

const T pipeLength = 7.4e-1; // in meters
const T globalInletPressure = 1.e07; // in Pascal or N/m^2
const T pressureGradient = globalInletPressure/pipeLength;
const T inletPressure = pressureGradient*lx0;

const T nu = 1.8800247371675945e-07;
const T nuLb = dx*dx/dt;

const T maxPhysT = 2e-8; // max. simulation time in s, SI unit
const T physInterval = 1e-10; // interval for the convergence check in s
const T residuum = 1e-5; // residuum for the convergence check

int nEvery = 100;
const T outputIter = maxPhysT/nEvery;
const T rampVTime = 1e3*dt;
const T tuner = 0.; //0.99; // for partialSlip only: 0->bounceBack, 1->freeSlip

T tau = 3*(nu/nuLb) + 0.5;

// Stores geometry information in form of material numbers
SuperGeometry2D<T> prepareGeometry( UnitConverter<T,DESCRIPTOR> const& converter )
{
OstreamManager clout( std::cout,”prepareGeometry” );
clout << “Prepare Geometry …” << std::endl;

// set number of cuboids/blocks
#ifdef PARALLEL_MODE_MPI
const int noOfCuboids = singleton::mpi().getSize();
#else
const int noOfCuboids = lx0/lx1;
#endif

// setup channel
Vector<T,2> extendChannel( lx0,ly0 );
Vector<T,2> originChannel( 0, 0 );
std::shared_ptr<IndicatorF2D<T>> channel = std::make_shared<IndicatorCuboid2D<T>>( extendChannel, originChannel );

// setup step
Vector<T,2> extendStep( lx0-lx1,ly1 );
// Vector<T,2> originStep( 0, 0 );
Vector<T,2> originStep( lx1, 0 );
std::shared_ptr<IndicatorF2D<T>> step = std::make_shared<IndicatorCuboid2D<T>>( extendStep, originStep );

CuboidGeometry2D<T>* cuboidGeometry = new CuboidGeometry2D<T>( *(channel-step), converter.getConversionFactorLength(), noOfCuboids );
HeuristicLoadBalancer<T>* loadBalancer = new HeuristicLoadBalancer<T>( *cuboidGeometry );
// Instantiation of a superGeometry
SuperGeometry2D<T> superGeometry( *cuboidGeometry, *loadBalancer, 2 );

// material numbers from zero to 2 inside geometry defined by indicator
superGeometry.rename( 0,2 );
superGeometry.rename( 2,1,1,1 );

Vector<T,2> extendBC( 0,ly0 );
Vector<T,2> originBC;
IndicatorCuboid2D<T> inflow( extendBC, originBC );
// Set material number for inflow

if(inletType == constantVel){
superGeometry.rename( 2,3,1,inflow );
}else if(inletType == constantPress){
superGeometry.rename( 2,6,1,inflow );
}else{}

originBC[0] = lx0;
IndicatorCuboid2D<T> outflow( extendBC, originBC );
// Set material number for outflow
superGeometry.rename( 2,4,1,outflow );

// Removes all not needed boundary voxels outside the surface
superGeometry.clean();

// Removes all not needed boundary voxels inside the surface
superGeometry.innerClean();
superGeometry.checkForErrors();
superGeometry.getStatistics().print();

clout << “Prepare Geometry … OK” << std::endl;
return superGeometry;
}

// Set up the geometry of the simulation
void prepareLattice( UnitConverter<T,DESCRIPTOR> const& converter,
SuperLattice2D<T,DESCRIPTOR>& sLattice,
Dynamics<T,DESCRIPTOR>& bulkDynamics,
sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& bc,
SuperGeometry2D<T>& superGeometry ) {
OstreamManager clout( std::cout,”prepareLattice” );
clout << “Prepare Lattice …” << std::endl;

const T omega = converter.getLatticeRelaxationFrequency();

auto bulkIndicator = superGeometry.getMaterialIndicator({1, 3, 4, 6});

// Material=0 –>do nothing
sLattice.defineDynamics( superGeometry, 0, &instances::getNoDynamics<T,DESCRIPTOR>() );
// Material=1 –>bulk dynamics
// Material=3 –>bulk dynamics (inflow)
// Material=4 –>bulk dynamics (outflow)
sLattice.defineDynamics( bulkIndicator, &bulkDynamics );

// Material=2 –>bounce-back || free-slip || partial-slip || no-slip
if (boundaryType == bounceBack) {
sLattice.defineDynamics( superGeometry, 2, &instances::getBounceBack<T, DESCRIPTOR>() );
} else if (boundaryType == freeSlip) {
sLattice.defineDynamics(superGeometry, 2, &instances::getNoDynamics<T, DESCRIPTOR>());
bc.addSlipBoundary( superGeometry, 2 );
// Partial Slip does not work with the superGeometry based creation of simulation domain,
// probably need to create boundary conditions and geometry explicitly.
} else if (boundaryType == partialSlip) {
sLattice.defineDynamics(superGeometry, 2, &instances::getNoDynamics<T, DESCRIPTOR>());
bc.addPartialSlipBoundary(tuner, superGeometry, 2 );
} else {
sLattice.defineDynamics( superGeometry, 2, &bulkDynamics );
bc.addVelocityBoundary( superGeometry, 2, omega );
}

// Initial conditions
AnalyticalConst2D<T,T> ux( 0. );
AnalyticalConst2D<T,T> uy( 0. );
AnalyticalConst2D<T,T> rho( 1. );
AnalyticalComposed2D<T,T> u( ux,uy );

//Initialize all values of distribution functions to their local equilibrium
sLattice.defineRhoU( bulkIndicator, rho, u );
sLattice.iniEquilibrium( bulkIndicator, rho, u );

// Setting of the boundary conditions
if(inletType == constantVel){
bc.addVelocityBoundary( superGeometry, 3, converter.getLatticeRelaxationFrequency() );
}else if(inletType == constantPress){
bc.addPressureBoundary( superGeometry, 6, converter.getLatticeRelaxationFrequency() );
AnalyticalConst2D<T,T> inFlowRhoBC(converter.getLatticeDensityFromPhysPressure( inletPressure ));
sLattice.defineRho( superGeometry, 6, inFlowRhoBC );
}else{}

bc.addPressureBoundary( superGeometry, 4, converter.getLatticeRelaxationFrequency() );
AnalyticalConst2D<T,T> outFlowRhoBC(converter.getLatticeDensityFromPhysPressure( outletPressure ));
sLattice.defineRho( superGeometry, 4, outFlowRhoBC );

// Make the lattice ready for simulation
sLattice.initialize();

clout << “Prepare Lattice … OK” << std::endl;
}

// Generates a slowly increasing inflow for the first iTMaxStart timesteps
void setBoundaryValues( UnitConverter<T,DESCRIPTOR> const& converter,
SuperLattice2D<T,DESCRIPTOR>& sLattice, int iT,
SuperGeometry2D<T>& superGeometry ) {
OstreamManager clout( std::cout,”setBoundaryValues” );

// clout << “I got here …” << std::endl;
// time for smooth start-up
int iTmaxStart = converter.getLatticeTime(rampVTime); //( maxPhysT*0.2 );
int iTupdate = 5;

if ( iT%iTupdate == 0 && iT<= iTmaxStart ) {
// clout << “I got alright here …” << std::endl;
// Smooth start curve, sinus
// SinusStartScale<T,int> StartScale(iTmaxStart, T(1));
// Smooth start curve, polynomial
PolynomialStartScale<T,int> StartScale( iTmaxStart, T( 1 ) );
// Creates and sets the Poiseuille inflow profile using functors
int iTvec[1] = {iT};
T frac[1] = {};
StartScale( frac,iTvec );
T maxVelocity = converter.getCharLatticeVelocity()*3./2.*frac[0];
T distance2Wall = converter.getConversionFactorLength()/2.;
Poiseuille2D<T> poiseuilleU( superGeometry, 3, maxVelocity, distance2Wall );
// define lattice speed on inflow
sLattice.defineU( superGeometry, 3, poiseuilleU );
}
}

// write data to termimal and file system
void getResults( SuperLattice2D<T,DESCRIPTOR>& sLattice,
UnitConverter<T,DESCRIPTOR> const& converter, int iT,
SuperGeometry2D<T>& superGeometry, Timer<T>& timer,
SuperPlaneIntegralFluxVelocity2D<T>& velocityFlux,
SuperPlaneIntegralFluxPressure2D<T>& pressureFlux, bool hasConverged ) {
OstreamManager clout( std::cout,”getResults” );
SuperVTMwriter2D<T> vtmWriter( “fstep2d” );

if ( iT==0 ) {
// Writes geometry, cuboid no. and rank no. to file system
SuperLatticeGeometry2D<T,DESCRIPTOR> geometry( sLattice, superGeometry );
SuperLatticeCuboid2D<T,DESCRIPTOR> cuboid( sLattice );
SuperLatticeRank2D<T,DESCRIPTOR> rank( sLattice );
vtmWriter.write( geometry );
vtmWriter.write( cuboid );
vtmWriter.write( rank );
vtmWriter.createMasterFile();
}

// Writes every outputIter seconds
if ( iT%converter.getLatticeTime( outputIter )==0 || hasConverged) {
SuperLatticePhysVelocity2D<T,DESCRIPTOR> velocity( sLattice, converter );
SuperLatticePhysPressure2D<T,DESCRIPTOR> pressure( sLattice, converter );
vtmWriter.addFunctor( velocity );
vtmWriter.addFunctor( pressure );
// write vtk to file system
vtmWriter.write( iT );
sLattice.communicate();
SuperEuklidNorm2D<T,DESCRIPTOR> normVel( velocity );
BlockReduction2D2D<T> planeReduction( normVel, 600, BlockDataSyncMode::ReduceOnly );
// write output as JPEG
heatmap::write(planeReduction, iT);
}

// Writes every outputIter*0.5 simulated
if ( iT%converter.getLatticeTime( outputIter*0.5 )==0 ) {
velocityFlux.print();
pressureFlux.print();

// write to terminal
timer.update( iT );
timer.printStep();
// Lattice statistics console output
sLattice.getStatistics().print( iT,converter.getPhysTime( iT ) );
}

// Saves lattice data
if ( iT%converter.getLatticeTime( 1 )==0 && iT>0 ) {
clout << “Checkpointing the system at t=” << iT << std::endl;
sLattice.save( “fstep2d.checkpoint” );
// The data can be reloaded using
// sLattice.load(“fstep2d.checkpoint”);
}
}

int main( int argc, char* argv[] ) {
// === 1st Step: Initialization ===
olbInit( &argc, &argv );
singleton::directories().setOutputDir( “./” ); // set output directory
OstreamManager clout( std::cout, “main” );

// create and open output file for saving input params
std::ofstream writer;
char fileName[40];
sprintf( fileName, “inputParams.dat”);
writer.open( fileName );
writer << std::setprecision( 10 );

writer << “All units are SI units.” << std::endl;
writer << ” boundaryType = ” << boundaryType << std::endl;
if(boundaryType == partialSlip){
writer << ” fraction of unchanged tangential flux along the channel walls, tuner = ” << tuner << std::endl;
}
writer << ” inletType = ” << inletType << std::endl;
writer << ” step position, lx1 = ” << lx1 << std::endl;
writer << ” step height, ly1 = ” << ly1 << std::endl;
writer << ” channel length, lx0 = ” << lx0 << std::endl;
writer << ” channel height, ly0 = ” << ly0 << std::endl;
writer << ” resolution of the model, N = ” << N << std::endl;
writer << ” lattice spacing, dx = ” << dx << std::endl;
writer << ” time step, dt = ” << dt << std::endl;
writer << ” fluid density, density = ” << density << std::endl;
if( inletType == constantVel ){
writer << ” inlet velocity BC, inletType = ” << uMax << std::endl;
}else if( inletType == timeVaryingVel ){
writer << ” ramp time for the time-varying inlet velocity BC, rampVTime = ” << rampVTime << std::endl;
}else if( inletType == constantPress ){
writer << ” Assumed nanofactory pipe Length = ” << pipeLength << std::endl;
writer << ” globalInletPressure = ” << globalInletPressure << std::endl;
writer << ” inlet pressure BC, inletPressure = ” << inletPressure << std::endl;
writer << ” total pressure gradient across the length of the channel, = ” << pressureGradient << std::endl;
}else{
clout << “Invalid –> inletType” << std::endl;
}
writer << ” outlet pressure BC, outletPressure = ” << outletPressure << std::endl;
writer << ” fluid viscosity, nu = ” << nu << std::endl;
writer << ” lattice fluid viscosity, nuLb = ” << nuLb << std::endl;
writer << ” lattice relaxation time, tau = ” << tau << std::endl;
writer << ” maximum simulation time, maxPhysT = ” << maxPhysT << std::endl;
writer << ” simulation time between writing outputs, outputIter = ” << outputIter << std::endl;
writer << ” simulation time between convergence checks, physInterval = ” << physInterval << std::endl;
writer << ” residual for convergence check, residuum = ” << residuum << std::endl;

writer.close();

// UnitConverter<T,DESCRIPTOR> converter(
// (T) 1./N, // physDeltaX: spacing between two lattice cells in __m__
// (T) 1./(M*N), // physDeltaT: time step in __s__
// (T) 1., // charPhysLength: reference length of simulation geometry
// (T) 0.5, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
// (T) 1./500., // physViscosity: physical kinematic viscosity in __m^2 / s__
// (T) 1. // physDensity: physical density in __kg / m^3__
// );

// if (argc > 1) {
// if (argv[1][0]==’-‘&&argv[1][1]==’h’) {
// OstreamManager clout( std::cout,”help” );
// clout<<“Usage: program [Resolution] [FlowType] [BoundaryType]”<<std::endl;
// clout<<“FlowType: 0=forced, 1=nonForced”<<std::endl;
// clout<<“BoundaryType: 0=bounceBack, 1=local, 2=interpolated, 3=freeSlip, 4=partialSlip”<<std::endl;
// clout<<“Default: FlowType=forced, Resolution=50, BoundaryType=interpolated”<<std::endl;
// return 0;
// }
// }

UnitConverterFromResolutionAndRelaxationTime<T, DESCRIPTOR> const converter(
int {N}, // resolution: number of voxels per charPhysL
(T) tau, // latticeRelaxationTime: relaxation time, have to be greater than 0.5!
(T) ly1, //1, // charPhysLength: reference length of simulation geometry
(T) uMax, //1, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) nu, //1./Re, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) density // physDensity: physical density in __kg / m^3__
);

// Prints the converter log as console output
converter.print();
// Writes the converter log in a file
converter.write(“fstep2d”);

// === 2nd Step: Prepare Geometry ===
// Instantiation of a superGeometry
SuperGeometry2D<T> superGeometry( prepareGeometry(converter) );

// === 3rd Step: Prepare Lattice ===
SuperLattice2D<T,DESCRIPTOR> sLattice( superGeometry );

BGKdynamics<T,DESCRIPTOR> bulkDynamics (
converter.getLatticeRelaxationFrequency(),
instances::getBulkMomenta<T,DESCRIPTOR>()
);

// Dynamics<T, DESCRIPTOR>* bulkDynamics;
//#if defined(MRT)
// if (flowType == forced) {
// ForcedMRTdynamics<T, DESCRIPTOR> bulkDynamics( converter.getLatticeRelaxationFrequency(), instances::getBulkMomenta<T, DESCRIPTOR>() );
// } else {
// bulkDynamics = new MRTdynamics<T, DESCRIPTOR> ( converter.getLatticeRelaxationFrequency(), instances::getBulkMomenta<T, DESCRIPTOR>() );
// }
//#else
// if (flowType == forced) {
// ForcedBGKdynamics<T, DESCRIPTOR> bulkDynamics( converter.getLatticeRelaxationFrequency(), instances::getBulkMomenta<T, DESCRIPTOR>() );
// } else {
// BGKdynamics<T, DESCRIPTOR> bulkDynamics( converter.getLatticeRelaxationFrequency(), instances::getBulkMomenta<T, DESCRIPTOR>() );
// }
//#endif

// choose between local and non-local boundary condition
sOnLatticeBoundaryCondition2D<T,DESCRIPTOR> sBoundaryCondition( sLattice );
createInterpBoundaryCondition2D<T,DESCRIPTOR>( sBoundaryCondition );
// createLocalBoundaryCondition2D<T,DESCRIPTOR>( sBoundaryCondition );

prepareLattice( converter, sLattice, bulkDynamics, sBoundaryCondition, superGeometry );

// instantiate reusable functors
SuperPlaneIntegralFluxVelocity2D<T> velocityFlux( sLattice, converter, superGeometry, {19.,1.}, {0.,1.} );
SuperPlaneIntegralFluxPressure2D<T> pressureFlux( sLattice, converter, superGeometry, {19.,1.}, {0.,1.} );

// === 4th Step: Main Loop with Timer ===
clout << “starting simulation…” << std::endl;
Timer<T> timer( converter.getLatticeTime( maxPhysT ), superGeometry.getStatistics().getNvoxel() );
util::ValueTracer<T> converge( converter.getLatticeTime( physInterval ), residuum );
timer.start();

if(inletType == constantVel){
T maxVelocity = converter.getCharLatticeVelocity()*3./2.;
T distance2Wall = converter.getConversionFactorLength()/2.;
Poiseuille2D<T> poiseuilleU( superGeometry, 3, maxVelocity, distance2Wall );
// define lattice speed on inflow
sLattice.defineU( superGeometry, 3, poiseuilleU );
}

for ( int iT = 0; iT < converter.getLatticeTime( maxPhysT ); ++iT ) {
if ( converge.hasConverged() ) {
clout << “Simulation converged.” << std::endl;
getResults( sLattice, converter, iT, superGeometry, timer, velocityFlux, pressureFlux, converge.hasConverged() );

break;
}

if(inletType == timeVaryingVel){
// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues( converter, sLattice, iT, superGeometry );
}

// === 6th Step: Collide and Stream Execution ===
sLattice.collideAndStream();

// === 7th Step: Computation and Output of the Results ===
getResults( sLattice, converter, iT, superGeometry, timer, velocityFlux, pressureFlux, converge.hasConverged() );
converge.takeValue( sLattice.getStatistics().getAverageEnergy(), true );
}

timer.stop();
timer.printSummary();
}

[Adrian]

Thanks for your detailed report!

We are aware of this problem and it will be fixed in our next release.
See the following patch for a solution:

`
Date: Fri, 5 Jun 2020 16:54:38 +0200
Subject: [PATCH] Fix bstep2d material number setup

_Implicit_ geometry modelling by approximating the step shape using the
cuboids themselves only works if we have enought cuboids to start with.

The material numbers were set correctly in e.g. olb 1.2. It probably got
lost somewhere along the introduction of managed functor arithmetic.
—
examples/laminar/bstep2d/bstep2d.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)

diff –git a/examples/laminar/bstep2d/bstep2d.cpp b/examples/laminar/bstep2d/bstep2d.cpp
index f6359164e..83faa3093 100644
— a/examples/laminar/bstep2d/bstep2d.cpp
+++ b/examples/laminar/bstep2d/bstep2d.cpp
@@ -85,7 +85,7 @@ SuperGeometry2D<T> prepareGeometry( UnitConverter<T,DESCRIPTOR> const& converter
SuperGeometry2D<T> superGeometry( *cuboidGeometry, *loadBalancer, 2 );

// material numbers from zero to 2 inside geometry defined by indicator
– superGeometry.rename( 0,2 );
+ superGeometry.rename( 0,2, *(channel-step) );
superGeometry.rename( 2,1,1,1 );

Vector<T,2> extendBC( 0,ly0 );
—
2.25.4
`

It works for larger cuboid numbers as the block decomposition can then approximate the step geometry as the unification of multiple cuboids. For single cuboid setups we have to use material numbers.

[vthk]

Hi Adrian,

Thanks very much for your response and for posting the work-around for the bug.
I found one of my own also using the IndicatorF2D() and IndicatorIdentity2D() approach shown below.

Thanks once again for taking the time to respond.

Warm regards,
User-Vthk.

int main( int argc, char* argv[] ) {
…

// === 2nd Step: Prepare Geometry ===
// Instantiation of a superGeometry
std::vector<T> extend = { lx0, ly0 };
std::vector<T> origin = { 0, 0 };
IndicatorCuboid2D<T> cuboid(extend,origin);

#ifdef PARALLEL_MODE_MPI
CuboidGeometry2D<T> cGeometry( cuboid, converter.getPhysDeltaX(), singleton::mpi().getSize() );
#else
CuboidGeometry2D<T> cGeometry( cuboid, converter.getPhysDeltaX() );
#endif

// Instantiation of loadbalancer
HeuristicLoadBalancer<T> loadBalancer( cGeometry );
loadBalancer.print();

SuperGeometry2D<T> superGeometry( cGeometry,loadBalancer );
prepareGeometry(superGeometry, converter );

…
}

void prepareGeometry( SuperGeometry2D<T>& superGeometry, UnitConverter<T,DESCRIPTOR> const& converter )
{
…
// setup channel
Vector<T,2> extendChannel( lx0,ly0 );
Vector<T,2> originChannel( 0, 0 );
std::shared_ptr<IndicatorF2D<T>> channel = std::make_shared<IndicatorCuboid2D<T>>( extendChannel,
originChannel );

// setup step
Vector<T,2> extendStep( lx0-lx1,ly1 );
Vector<T,2> originStep( lx1, 0 );
std::shared_ptr<IndicatorF2D<T>> step = std::make_shared<IndicatorCuboid2D<T>>( extendStep,
originStep );

// setup boolean operation for the geometry
IndicatorIdentity2D<T> fstep2d( channel – step ); //

// material numbers from zero to 2 inside geometry defined by indicator
superGeometry.rename( 0,2,fstep2d );
superGeometry.rename( 2,1,1,1 );
…
}

[Adrian]

Happy to help.

As a sidenote: Your solution is basically the same as what is done in my patch – *(channel-step) is just a shorthand for wrapping the subtracted indicators in an identity.

[vthk]

Yes, in effect it is! 🙂
:thumbsup:

Regards,
User-Vthk.

[Author]

Posts