Importing stl in FreesurfaceBreakingDam model.

[Dip]

Hi,
I tried to modified the breakingDam3d, I want to input a stl in the already given geometry. I want to insert the geometry at the middle. But, my stl is not loading. What is the problem with it?

This is my code:
/* This file is part of the OpenLB library
*
* Copyright (C) 2021 Claudius Holeksa, Maximilian Schecher
* E-mail contact: info@openlb.net
* The most recent release of OpenLB can be downloaded at
* <http://www.openlb.net
*
* This program 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 2
* of the License, or (at your option) any later version.
*
* This program 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 this program; if not, write to the Free
* Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301, USA.
*/

#include “olb3D.h”
#include “olb3D.hh”

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

using namespace olb;
using namespace olb::descriptors;

using T = float;
using DESCRIPTOR = D3Q27<descriptors::FORCE, FreeSurface::MASS, FreeSurface::EPSILON, FreeSurface::CELL_TYPE, FreeSurface::CELL_FLAGS, FreeSurface::TEMP_MASS_EXCHANGE, FreeSurface::PREVIOUS_VELOCITY>;

/*
* Helper since there are a lot of values to set and giving a reference to this object is easier than
* defining the function calls accordingly
*/
struct FreeSurfaceAppHelper {
std::array<T,3> area{{1,0.5,0.5}};
std::array<T,3> gravity_force = {{0.,0., -9.81}};

T char_phys_length = 10;
T char_phys_vel = 0.1;
bool has_surface_tension = true;
T surface_tension_coefficient = 0.0661;

};

template<typename T, typename DESCRIPTOR>
class FreeSurfaceBreakingDam3D final : public AnalyticalF3D<T,T> {
private:
T lattice_size;
// This value doesn’t depend on the dimension of the descriptor. It’s always 3
std::array<T,4> cell_values;
std::array<T,3> area;

// STLreader member to check if a point is inside the STL
STLreader<T> stlReader;

public:
FreeSurfaceBreakingDam3D(T lattice_size_, const std::array<T,4>& cell_vals, const std::array<T,3>& area_)
: AnalyticalF3D<T,T>{1}, lattice_size{lattice_size_}, cell_values{cell_vals}, area{area_},
stlReader(“KCS_hull_SVA_0.001.STL”, 1e-3, 1, 2, false, 0., 1.) {}

// This sets which parts of the area are going to be the fluid
bool operator()(T output[], const T x[]) override {
output[0] = cell_values[0];

// Array to store whether the point is inside the STL
bool inside[1];

// Check if the point x is inside the STL geometry
if (stlReader(inside, x) && inside[0]) {
output[0] = cell_values[4]; // Set the cell value to 4 if inside the STL
} else if (x[2] <= area[2] * 0.5) {
output[0] = cell_values[2];
} else if (x[2] – lattice_size * 1.1 <= area[2] * 0.5 ) {
output[0] = cell_values[1];
}

return true;
}
};

void prepareGeometry( UnitConverter<T,DESCRIPTOR> const& converter,
SuperGeometry<T,3>& superGeometry, STLreader<T>& stlReader ) {

OstreamManager clout( std::cout, “prepareGeometry” );
clout << “Prepare Geometry …” << std::endl;

// Initial renaming of material numbers
superGeometry.rename( 0, 2 ); // Rename material 0 to 2
superGeometry.rename( 2, 1, {1,1,1} ); // Rename material 2 to 1 for fluid areas

// Rename cells inside the STL geometry to material 4
superGeometry.rename( 1, 4, stlReader );

// Cleaning and finalizing the geometry
superGeometry.clean();
superGeometry.innerClean();
superGeometry.checkForErrors();
superGeometry.print();

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

void prepareBreakingDam(UnitConverter<T,DESCRIPTOR> const& converter,
SuperLattice<T, DESCRIPTOR>& sLattice,
SuperGeometry<T,3>& superGeometry, T lattice_size, const FreeSurfaceAppHelper& helper)
{
AnalyticalConst3D<T,T> zero( 0. );
AnalyticalConst3D<T,T> one( 1. );
AnalyticalConst3D<T,T> two( 2. );
AnalyticalConst3D<T,T> four( 4. );
FreeSurfaceBreakingDam3D<T,DESCRIPTOR> cells_analytical{ lattice_size, {0., 1., 2., 4.}, helper.area};
FreeSurfaceBreakingDam3D<T,DESCRIPTOR> mass_analytical{ lattice_size, {0., 0.5, 1., 2.}, helper.area};

AnalyticalConst3D<T,T> force_zero{0., 0., 0.};

for (int i: {0,1,2,4}) {
sLattice.defineField<FreeSurface::MASS>(superGeometry, i, zero);
sLattice.defineField<FreeSurface::EPSILON>(superGeometry, i, zero);
sLattice.defineField<FreeSurface::CELL_TYPE>(superGeometry, i, zero);
sLattice.defineField<FreeSurface::CELL_FLAGS>(superGeometry, i, zero);
sLattice.defineField<descriptors::FORCE>(superGeometry, i, force_zero);
}

sLattice.defineField<FreeSurface::CELL_TYPE>(superGeometry, 1, cells_analytical);

sLattice.defineField<FreeSurface::MASS>(superGeometry, 1, mass_analytical);
sLattice.defineField<FreeSurface::EPSILON>(superGeometry, 1, mass_analytical);

for (int i: {0,2}) {
//sLattice.defineField<FreeSurface::MASS>(superGeometry, i, one);
sLattice.defineField<FreeSurface::EPSILON>(superGeometry, i, one);
sLattice.defineField<FreeSurface::CELL_TYPE>(superGeometry, i, four);
}

T force_factor = 1./ converter.getConversionFactorForce() * converter.getConversionFactorMass();
AnalyticalConst3D<T,T> force_a{helper.gravity_force[0] * force_factor, helper.gravity_force[1] * force_factor, helper.gravity_force[2] * force_factor};
sLattice.defineField<descriptors::FORCE>(superGeometry.getMaterialIndicator(1), force_a);

}

void prepareLattice( UnitConverter<T,DESCRIPTOR> const& converter,
SuperLattice<T, DESCRIPTOR>& sLattice,
SuperGeometry<T,3>& superGeometry, T lattice_size, const FreeSurfaceAppHelper& helper) {

OstreamManager clout( std::cout,”prepareLattice” );
clout << “Prepare Lattice …” << std::endl;

// Material=0 –>do nothing
sLattice.defineDynamics<NoDynamics<T,DESCRIPTOR>>(superGeometry, 0);
// Material=1 –>bulk dynamics
sLattice.defineDynamics<SmagorinskyForcedBGKdynamics<T,DESCRIPTOR>>( superGeometry, 1);
// Material=2 –>no-slip boundary
sLattice.defineDynamics<BounceBack<T,DESCRIPTOR>>( superGeometry, 2);
setSlipBoundary<T,DESCRIPTOR>(sLattice, superGeometry, 4);

sLattice.setParameter<descriptors::OMEGA>(converter.getLatticeRelaxationFrequency());
sLattice.setParameter<collision::LES::Smagorinsky>(T(0.2));

prepareBreakingDam(converter, sLattice, superGeometry, lattice_size, helper);

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

void setInitialValues(SuperLattice<T, DESCRIPTOR>& sLattice, SuperGeometry<T,3>& sGeometry, T lattice_length, UnitConverter<T,DESCRIPTOR> const& converter){
OstreamManager clout( std::cout,”setInitialValues” );

AnalyticalConst3D<T,T> u{0., 0.};
AnalyticalConst3D<T,T> one(1.);

sLattice.defineRhoU( sGeometry.getMaterialIndicator({0,1,2,4}), one, u );
for (int i: {0,1,2,4}) {
sLattice.iniEquilibrium( sGeometry, i, one, u );
}

// Set up free surface communicator stages
FreeSurface::initialize(sLattice);
// Make the lattice ready for simulation
sLattice.initialize();
}

void getResults( SuperLattice<T,DESCRIPTOR>& sLattice,
UnitConverter<T,DESCRIPTOR> const& converter, int iT,
SuperGeometry<T,3>& superGeometry, util::Timer<T>& timer)
{
const int vtmIter = 100;
const int statIter = 100;

if ( iT==0 ) {
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperVTMwriter3D<T> vtmWriter( “breakingDam3d” );
SuperLatticeGeometry3D<T, DESCRIPTOR> geometry( sLattice, superGeometry );
SuperLatticeCuboid3D<T, DESCRIPTOR> cuboid( sLattice );
SuperLatticeRank3D<T, DESCRIPTOR> rank( sLattice );
vtmWriter.write( geometry );
vtmWriter.write( cuboid );
vtmWriter.write( rank );

vtmWriter.createMasterFile();
}

// Writes output in ParaView
if ( iT%vtmIter==0 ) {
sLattice.setProcessingContext(ProcessingContext::Evaluation);

SuperVTMwriter3D<T> vtmWriter( “breakingDam3d” );
SuperLatticePhysVelocity3D<T, DESCRIPTOR> velocity( sLattice, converter );
SuperLatticePhysPressure3D<T, DESCRIPTOR> pressure( sLattice, converter );
SuperLatticeExternalScalarField3D<T, DESCRIPTOR, FreeSurface::EPSILON> epsilon( sLattice );
SuperLatticeExternalScalarField3D<T, DESCRIPTOR, FreeSurface::CELL_TYPE> cells( sLattice );
SuperLatticeExternalScalarField3D<T, DESCRIPTOR, FreeSurface::MASS> mass( sLattice );
SuperLatticeGeometry3D<T, DESCRIPTOR> geometry(sLattice, superGeometry);
epsilon.getName() = “epsilon”;
cells.getName() = “cell_type”;
mass.getName() = “mass”;
vtmWriter.addFunctor( velocity );
vtmWriter.addFunctor( pressure );
vtmWriter.addFunctor( epsilon );
vtmWriter.addFunctor( cells );
vtmWriter.addFunctor( mass );
vtmWriter.addFunctor( geometry );

vtmWriter.write(iT);
}

// Writes output on the console
if ( iT%statIter==0 ) {
// Timer console output
timer.update( iT );
timer.printStep();

// Lattice statistics console output
sLattice.getStatistics().print( iT,converter.getPhysTime( iT ) );
}
}

namespace {
FreeSurfaceAppHelper free_surface_config {
{10.,.5,.5}, //Area
{0.,0.,-9.81}, // Gravity force
10., // char_phys_length
0.1, // char_phys_vel
true, // has_surface_tension
0.05 // surface tension coefficient
};

class FreeSurfaceConfig {
public:
T viscosity = 1e-4; // viscosity of water: 1e-4
//T viscosity = 122000e-4; // viscosity of honey: 122000e-4 = 12.2
T density = 1e3; // density of water: 1e3
//T density = 1.358e3; // density of honey: 1.358e3 = 1358
T physTime = 40.;
T latticeRelaxationTime = .501;
int N = 500;

// Anti jitter value
T transitionThreshold = 1e-3;
// When to remove lonely cells
T lonelyThreshold = 1.0;
};

}

int main(int argc, char **argv)
{
olbInit(&argc, &argv, false, false);

FreeSurfaceConfig c;
OstreamManager clerr( std::cerr, “main” );
OstreamManager clout( std::cout, “main” );

singleton::directories().setOutputDir(“./tmp/”);

FreeSurfaceAppHelper& helper = free_surface_config;

UnitConverterFromResolutionAndRelaxationTime<T, DESCRIPTOR> const converter(
int {c.N}, // resolution: number of voxels per charPhysL
(T) c.latticeRelaxationTime, // latticeRelaxationTime: relaxation time, have to be greater than 0.5!
(T) helper.char_phys_length, // charPhysLength: reference length of simulation geometry
(T) helper.char_phys_vel, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) c.viscosity, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) c.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(“free surface”);

T lattice_size = helper.char_phys_length / c.N;

T force_conversion_factor = 1./converter.getConversionFactorForce()*converter.getConversionFactorMass();

// Convert kg / s^2
// Basically it is multiplied with s^2 / kg = s^2 * m^3 / (kg * m^2 * m) = 1. / (velocity_factor^2 * density * length_factor)
T surface_tension_coefficient_factor = std::pow(converter.getConversionFactorTime(),2)/ (c.density * std::pow(converter.getConversionFactorLength(),3));

clout<<“Surface: “<<surface_tension_coefficient_factor * helper.surface_tension_coefficient<<std::endl;
clout<<“Lattice Size: “<<converter.getPhysDeltaX()<<std::endl;

// === 2nd Step: Prepare Geometry ===
Vector<T,3> extend( helper.area[0], helper.area[1], helper.area[2] );
Vector<T,3> origin;
IndicatorCuboid3D<T> cuboid( extend, origin );

// Instantiation of a cuboidGeometry with weights
#ifdef PARALLEL_MODE_MPI
const int noOfCuboids = singleton::mpi().getSize();
#else
const int noOfCuboids = 4;
#endif
CuboidGeometry3D<T> cuboidGeometry( cuboid, converter.getConversionFactorLength(), noOfCuboids );

HeuristicLoadBalancer<T> loadBalancer( cuboidGeometry );
SuperGeometry<T,3> superGeometry( cuboidGeometry, loadBalancer, 2 );

STLreader<T> stlReader(“KCS_hull_SVA_0.001.STL”, 1e-3, 1, 2, false, 0., 1.);

prepareGeometry( converter, superGeometry,stlReader );

SuperLattice<T, DESCRIPTOR> sLattice( superGeometry );

clout<<“Overlap: “<<sLattice.getOverlap()<<std::endl;

prepareLattice( converter, sLattice, superGeometry, lattice_size, helper);

FreeSurface3DSetup<T,DESCRIPTOR> free_surface_setup{sLattice};

free_surface_setup.addPostProcessor();

// Set variables from freeSurfaceHelpers.h
sLattice.setParameter<FreeSurface::DROP_ISOLATED_CELLS>(true);
sLattice.setParameter<FreeSurface::TRANSITION>(c.transitionThreshold);
sLattice.setParameter<FreeSurface::LONELY_THRESHOLD>(c.lonelyThreshold);
sLattice.setParameter<FreeSurface::HAS_SURFACE_TENSION>(helper.has_surface_tension);
sLattice.setParameter<FreeSurface::SURFACE_TENSION_PARAMETER>(surface_tension_coefficient_factor * helper.surface_tension_coefficient);
sLattice.setParameter<FreeSurface::FORCE_CONVERSION_FACTOR>(force_conversion_factor);
sLattice.setParameter<FreeSurface::LATTICE_SIZE>(converter.getPhysDeltaX());

// === 4th Step: Main Loop with Timer ===
clout << “starting simulation…” << std::endl;
util::Timer<T> timer( converter.getLatticeTime( c.physTime ), superGeometry.getStatistics().getNvoxel() );
timer.start();
setInitialValues(sLattice, superGeometry, lattice_size, converter);

for ( std::size_t iT = 0; iT < converter.getLatticeTime( c.physTime ); ++iT ) {
getResults( sLattice, converter, iT, superGeometry, timer );
sLattice.collideAndStream();
}

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

return 0;
}

[Adrian]

What do you mean by “is not loading”?

One common issue when importing STLs is the scaling between “STL units” and the physical lengths used by the geometry setup in OpenLB.

[Dip]

Hello,
Thanks for the suggestion, but I tried scaling the stl while input 1, 1e-3,1e-6. But, after simulation results are like there’s no stl. How to fix this.

[Adrian]

Did you verify that the STL is correctly voxelized in the material geometry (i.e. that there are material 4 cells (going by your listing) where the STL should be)?

[Dip]

Sorry, I didn’t get what you mean by verifying. The stl is not developed in my setup model, when viewed in paraview. I mean there is no geometry of my stl.

[Adrian]

You can visualize the geometry in Paraview (I don’t mean the simulation results, I mean the geometry file that is also written out into the tmp folder). Another indication is if cells with the material number assigned to the STL geometry are listed in the geometry statistics that are printed out into the terminal when you run your case.

[Dip]

No, the stl is not loaded in the geometry.

[Adrian]

Ok, so you should take a closer look at your geometry setup in the prepareGeometry function.

[Dip]

Thanks. That’s my question. Any suggestion what might be wrong in my prepareGeometry function?

[Adrian]

No, without the geometry and running the case on my own (which I won’t do – you have to do this work yourself) I can only guess. With the given information I suspect that the scaling and/or placement of the STL are the problem.

[Dip]

Ok, thanks for the suggestions.

[Author]

Posts