Hello,

Here’s my function which defines the geometry!

Sorry it’s long but it’s essentially geomtry definition.

Any help would be greatly appreciated.

Mathieu

const T x1 = 0.;

const T y1l = -4.11e-4;

const T y1u = 4.11e-4;

const T x2 = 1.41e-3;

const T y2l = -4.11e-4;

const T y2u = 4.11e-4;

const T x3 = 4.e-3;

const T y3l = -2.6e-3;

const T y3u = 2.6e-3;

const T x4 = 2e-2;

const T y4l = -4.11e-4;

const T y4u = 4.11e-4;

const T L = 1.91e-2;

const T L0 = 8.22e-4; // hydraulic diameter

T poiseuilleVelocity(int iY, int N, T u) {

T y = (T)iY / (T)N;

return 4.*u * (y-y*y);

}

void iniGeometry( BlockStructure2D<T,DESCRIPTOR>& lattice,

LBunits const& converter,

Dynamics<T, DESCRIPTOR>& bulkDynamics,

OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition )

{

const T omega = converter.getOmega();

```
const T x0 = 0.;
const T y0 = -3.2e-3;
const T y_or = -y0/L0;
const T x1_s = (x1-x0)/L0;
const T y1l_s = (y1l-y0)/L0;
const T y1u_s = (y1u-y0)/L0;
const T x2_s = (x2-x0)/L0;
const T y2l_s = (y2l-y0)/L0;
const T y2u_s = (y2u-y0)/L0;
const T x3_s = (x3-x0)/L0;
const T y3l_s = (y3l-y0)/L0;
const T y3u_s = (y3u-y0)/L0;
const T x4_s = (x4-x0)/L0;
const T y4l_s = (y4l-y0)/L0;
const T y4u_s = (y4u-y0)/L0;
const T L_s = L/L0;
const int nx = converter.getNx();
const int ny = converter.getNy();
const int n_x0 = converter.nCell(x0);
const int n_y0 = converter.nCell(y_or);
const int n_x1 = converter.nCell(x1_s);
const int n_y1l = converter.nCell(y1l_s);
const int n_y1u = converter.nCell(y1u_s);
const int n_x2 = converter.nCell(x2_s);
const int n_y2l = converter.nCell(y2l_s);
const int n_y2u = converter.nCell(y2u_s);
const int n_x3 = converter.nCell(x3_s);
const int n_y3l = converter.nCell(y3l_s);
const int n_y3u = converter.nCell(y3u_s);
const int n_x4 = converter.nCell(x4_s);
const int n_y4l = converter.nCell(y4l_s);
const int n_y4u = converter.nCell(y4u_s);
const int n_L = converter.nCell(L_s);
T *js, *jn;
js = new T[nx];
jn = new T[nx];
// Bulk dynamics
lattice.defineDynamics(0,nx-1,0,ny-1, &bulkDynamics);
// Definition of the obstacle (bounce-back nodes)
int iY, nYmin, nYmax, iterate, iswitch=0, isw, inw, ise, ine;
T u[2] = {converter.getLatticeU(),0.};
T rho = (T)1;
cout << "n_y4l = " << n_y4l << "\n";
cout << "n_y4u = " << n_y4u << "\n";
cout << "n_y0 = " << n_y0 << "\n";
cout << "n_L = " << n_L << "\n";
cout << "n_x4 = " << n_x4 << "\n";
// Find the indices of upper and lower nodes on the wall
for (int iX=0; iX<nx; ++iX) {
iY = 0;
iterate = 1;
while (iterate==1)
{
if (iX<n_x2)
{
nYmin = n_y1l;
nYmax = n_y1u;
if ((iY>nYmin)&&(iY<nYmax))
{
if (iswitch==0)
{
iswitch=1;
js[iX] = iY;
if (iX==0)
{
isw = iY;
}
}
}
else
{
if (iswitch==1)
{
iswitch=0;
jn[iX] = iY-1;
if (iX==0)
{
inw = iY-1;
}
}
}
}
else if ((iX>=n_x2)&&(iX<n_x3))
{
nYmin = (n_y3l-n_y2l)*(iX-n_x2)/(n_x3-n_x2)+n_y2l;
nYmax = (n_y3u-n_y2u)*(iX-n_x2)/(n_x3-n_x2)+n_y2u;
if ((iY>nYmin)&&(iY<nYmax))
{
if (iswitch==0)
{
iswitch=1;
js[iX] = iY;
}
}
else
{
if (iswitch==1)
{
iswitch=0;
jn[iX] = iY-1;
}
}
}
else
{
nYmin = (n_y4l-n_y0)*n_L/(n_L-(n_x4-iX))+n_y0;
nYmax = (n_y4u-n_y0)*n_L/(n_L-(n_x4-iX))+n_y0;
if ((iY>nYmin)&&(iY<nYmax))
{
if (iswitch==0)
{
iswitch=1;
js[iX] = iY;
if (iX==(nx-1))
{
ise = iY;
}
}
}
else
{
if (iswitch==1)
{
iswitch=0;
jn[iX] = iY-1;
if (iX==(nx-1))
{
ine = iY-1;
}
}
}
}
if (iY>=(ny-1))
iterate = 0;
else
++iY;
}
}
// Define parabolic profile in the fluid
// Define bounce-back at the boundary
for (int iX=0; iX<nx; ++iX) {
iY = 0;
iterate = 1;
while (iterate==1)
{
if ((iY>=js[iX])&&(iY<=jn[iX]))
{
T vel[] = {
poiseuilleVelocity(iY-js[iX], jn[iX]-js[iX], converter.getLatticeU()),
T()
};
lattice.get(iX,iY).defineRhoU((T)1, vel);
lattice.get(iX,iY).iniEquilibrium((T)1, vel);;
}
else
{
lattice.defineDynamics( iX,iX,iY,iY,
&instances::getBounceBack<T,DESCRIPTOR>() );
}
if (iY>=(ny-1))
iterate = 0;
else
++iY;
}
}
cout << "isw = " << isw << "\n";
cout << "inw = " << inw << "\n";
cout << "ise = " << ise << "\n";
cout << "ine = " << ine << "\n";
boundaryCondition.addVelocityBoundary0N(0,0, isw,inw, omega);
boundaryCondition.addVelocityBoundary0P(nx-1,nx-1, ise, ine, omega);
// Make the lattice ready for simulation
lattice.initialize();
delete [] js;
delete [] jn;
```

}