Skip to content

Incorrect geometry creation – forward-facing 2D step

OpenLB – Open Source Lattice Boltzmann Code Forums on OpenLB Bug Reports Incorrect geometry creation – forward-facing 2D step

Viewing 5 posts - 1 through 5 (of 5 total)
  • Author
    Posts
  • #5130
    vthk
    Participant

    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();
    }

    #5131
    Adrian
    Participant

    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.

    #5133
    vthk
    Participant

    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 );

    }

    #5134
    Adrian
    Participant

    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.

    #5135
    vthk
    Participant

    Yes, in effect it is! 🙂
    :thumbsup:

    Regards,
    User-Vthk.

Viewing 5 posts - 1 through 5 (of 5 total)
  • You must be logged in to reply to this topic.