How to create 2D triangle ？

[Fany]

Hi everyone,

I want to create some 2D triangles but there is no IndicatorF for triangles like IndicatorCuboid2D for cubiods. Do you know how to do that?

I found there is a functor called SmoothIndicatorTriangle2D but what do the mass, epsilon and theta in that functor mean? The “SmoothXXX” one seems to be different from the general IndicatorF like IndicatorCuboid2D. What does “Smooth” stand for? Thank you.

[Adrian]

Does IndicatorTriangle2D fullfill your requirements?

This will be available in the next release. Ordinarily I’d quickly backport it to 1.4 for you but the new indicator is tied to more comprehensive changes of the library and unfortunately not easy to extract.

In general, smooth indicators are used by the resolved particle code to model particles with a porous boundary layer (in HLBM). See e.g.

R. Trunk et al. “Towards the Simulation of arbitrarily shaped 3D particles using a homogenised lattice Boltzmann method”. In: Computers & Fluids 172 (2018). doi: https://doi.org/10.1016/j.compfluid.2018.02.027.

[Fany]

Hi Adrian,

Glad to hear the new release. By the way, in the present vision 1.4, Does the parallelepiped version of creating a 2D cuboid achieve by setting the parameter Theta? If does, the the code below should include the parameter theta?

template <typename S>
IndicatorCuboid2D<S>::IndicatorCuboid2D(S xLength, S yLength, Vector<S,2> center, S theta )
  : _center(center), _xLength(xLength), _yLength(yLength), _theta(-theta)
{
  this->_myMin = {_center[0] - _xLength/2., _center[1] - _yLength/2.};
  this->_myMax = {_center[0] + _xLength/2., _center[1] + _yLength/2.};
}

[Fany]

Hi Adrian, I am creating the 2D triangle and cuboid without parallel to the X=0 and Y=0. But the computation with my geometry diverged. Could you please help check what the problems with the code below are?

// creating the cuboid with the theta of rotation around its center
template <typename S>
IndicatorCuboid2D<S>::IndicatorCuboid2D(S xLength, S yLength, Vector<S,2> center, S theta )
  : _center(center), _xLength(xLength), _yLength(yLength), _theta(theta)
{
  if (_theta>=0 && _theta <= M_PI2){
     this->_myMin = {_center[0] - (_xLength*std::cos(_theta)/2.+_yLength*std::sin(_theta)/2.), _center[1] - (_yLength*std::cos(_theta)/2.+_xLength*std::sin(_theta)/2.) };
     this->_myMax = {_center[0] + (_xLength*std::cos(_theta)/2.+_yLength*std::sin(_theta)/2.), _center[1] + (_yLength*std::cos(_theta)/2.+_xLength*std::sin(_theta)/2.) };
  } 
  if (_theta>=-M_PI2 && _theta < 0){
	 this->_myMin = {_center[0] - (_xLength*std::cos(_theta)/2.-_yLength*std::sin(_theta)/2.), _center[1] - (_yLength*std::cos(_theta)/2.-_xLength*std::sin(_theta)/2.) };
     this->_myMax = {_center[0] + (_xLength*std::cos(_theta)/2.-_yLength*std::sin(_theta)/2.), _center[1] + (_yLength*std::cos(_theta)/2.-_xLength*std::sin(_theta)/2.)};  
  }
}
template <typename S>
bool IndicatorCuboid2D<S>::operator()(bool output[], const S input[])
{
  S x, y;
  if ( !util::nearZero(_theta) ) {
    x = _center[0] + (input[0] - _center[0])*std::cos(_theta) + (input[1] - _center[1])*std::sin(_theta);  
    y = _center[1] + (input[1] - _center[1])*std::cos(_theta) - (input[0] - _center[0])*std::sin(_theta);
  } else {
    x = input[0];
    y = input[1];
  }

  output[0] = (  (fabs(_center[0] - x) < _xLength/2. || util::approxEqual(fabs(_center[0] - x),_xLength/2.) )
                 && (fabs(_center[1] - y) < _yLength/2. || util::approxEqual(fabs(_center[1] - y), _yLength/2.) ) );
  return true;
}

//  creating the triangle
template <typename S>
IndicatorTriangle2D<S>::IndicatorTriangle2D(Vector<S,2> PointA, Vector<S,2> PointB, Vector<S,2> PointC)
 : _PointA(PointA), _PointB(PointB), _PointC(PointC)
{
   
  this->_myMin = { std::min(_PointA[0], std::min(_PointB[0], _PointC[0])),  std::min(_PointA[1], std::min(_PointB[1], _PointC[1]))};
  this->_myMax = { std::max(_PointA[0], std::max(_PointB[0], _PointC[0])),  std::max(_PointA[1], std::max(_PointB[1], _PointC[1]))};

  _ab = _PointB - _PointA;
  _bc = _PointC - _PointB;
  _ca = _PointA - _PointC;
  
}

template <typename S>
bool IndicatorTriangle2D<S>::operator()(bool output[], const S input[])
{
  Vector<S, 2> _Point ;
  _Point[0] = input[0];
  _Point[1] = input[1]; //any point M
  Vector<S, 2> _am=_Point - _PointA ; 
  Vector<S, 2> _bm=_Point - _PointB ;
  Vector<S, 2> _cm=_Point - _PointC ;
    
  S S_ABM = fabs(_ab[1]*_am[0]-_ab[0]*_am[1])/2. ;   // the area of triangle ABM;
  S S_BCM = fabs(_bc[1]*_bm[0]-_ca[0]*_bm[1])/2. ;  // the area of triangle BCM;
  S S_CAM = fabs(_ca[1]*_cm[0]-_ca[0]*_cm[1])/2. ;  // the area of triangle CAM;
  S Sum = fabs(_ca[1]*_ab[0]-_ca[0]*_ab[1])/2. ;  // the area of triangle CAM;
  
  output[0] = util::approxEqual(S_ABM + S_BCM + S_CAM,Sum) ;  
  return true;   
}

[Adrian]

So the simulation worked with a different indicator? Do the material numbers look as expected (you can check in ParaView)? Which boundary conditions and bulk dynamics are you using? If you inspect the results prior to divergence, where do the problems originate from?

[Fany]

Yes, I created a new indicator for triangle and inclined rectangle, but i don’t know why there is a problem in setting the outlet materials for inclined rectangle and how to clean the unnecessary voxels outside, as following:

The statistic of outflow material number (another geometry of inclined rectangle) was 1-fluid, 2-boundary, 3-inflow, 4-outflow:
`[SuperGeometryStatistics2D]materialNumber=1; count=18848; minPhysR=(-0.00148286,-0.00185643); maxPhysR=(0.00324857,0.00186357)
[SuperGeometryStatistics2D] materialNumber=2; count=1066; minPhysR=(-0.0015,-0.00187357); maxPhysR=(0.00326571,0.00188072)
[SuperGeometryStatistics2D] materialNumber=3; count=45; minPhysR=(-0.0015,-0.000382141); maxPhysR=(-0.0015,0.000372145)
[SuperGeometryStatistics2D] materialNumber=4; count=3; minPhysR=(0.00321429,0.00148643); maxPhysR=(0.00326571,0.00153786)

[Adrian]

This looks as if the inclined rectangle intersects the domain boundaries which likely interferes with the setup. I would also suggest to increase the resolution as I would be surprised if you got sensible results for this rough a setup.

[Fany]

It is not resolution but the number of cuboids ( by CuboidGeometry). You can see the following image.

[Adrian]

Ok – in this case there is definitely something weird going on as the rectangle intersects the cuboids and these intersections are also where the divergence seems to come from. What happens if you use only one cuboid (as a minimal case)? (In any case: For such a low resolution it likely doesn’t make sense to increase the cuboid number to the degree that it is done here)

Not directly related: What is the reason for inclining the rectangle if you only simulate inside it anyway?

• This reply was modified 2 years, 5 months ago by Adrian.

[Fany]

Good, thanks. I will try to increase the resolution. But I found the resolution in the examples of OpenLB library seemed to be usually set at 10-35. In my case, the resolution was set at 10 (same with that in rayleiBenard2d case).
By the way, the inclined rectangle above was just obtained by rotating the rectangle of rayleiBenard2d case (OpenLB) for 45°. The blue area above should be the unnecessary part outside of the rectangle?

[Fany]

I am trying to simulate the flow heat transfer in a inclined pipe (2D)

[Adrian]

w.r.t. resolution: The N constant (which as you say is commonly between 10 – 35) is only a useful proxy for the actual resolution that results of the spatial discretization delta and the domain size. This delta is computed by dividing some reference length (chosen to fit what is modeled by each simulation case) divided by the N constant.

w.r.t. inclined geometry / pipe: If my understanding of what you want to simulate is correct, you would be better served to keep the geometry axis aligned s.t. it is nicely discretized by the regular grid and model the inclination by rotating the external force field. This also should get rid of the cuboid decomposition problems (which point to some issue in domain setup).

[Fany]

Hi Adrian, I was creating the 2D triangle by my own triangle functor. But the visualization result showed my functor did not work. What could be the problems here? By the way, when will the next release be ？
`
// creating the triangle
template <typename S>
IndicatorTriangle2D<S>::IndicatorTriangle2D(Vector<S,2> PointA, Vector<S,2> PointB, Vector<S,2> PointC)
: _PointA(PointA), _PointB(PointB), _PointC(PointC)
{
// S AB=std::sqrt(std::pow(_PointB[0]-_PointA[0], 2)+std::pow(_PointB[1]-_PointA[1], 2));
// S AC=std::sqrt(std::pow(_PointC[0]-_PointA[0], 2)+std::pow(_PointC[1]-_PointA[1], 2));
// S BC=std::sqrt(std::pow(_PointB[0]-_PointC[0], 2)+std::pow(_PointB[1]-_PointC[1], 2));
// S halfperi=0.5*(AB+AC+BC);
// S temp=halfperi*(halfperi-AB)*(halfperi-AC)*(halfperi-BC);
// this->_circumRadius = 0.25*AB*AC*BC/std::sqrt(temp);

this->_myMin = { std::min(_PointA[0], std::min(_PointB[0], _PointC[0])), std::min(_PointA[1], std::min(_PointB[1], _PointC[1]))};
this->_myMax = { std::max(_PointA[0], std::max(_PointB[0], _PointC[0])), std::max(_PointA[1], std::max(_PointB[1], _PointC[1]))};
// this->_theta = theta * M_PI/180.;

_ab = _PointB – _PointA;
_bc = _PointC – _PointB;
_ca = _PointA – _PointC;

}

template <typename S>
bool IndicatorTriangle2D<S>::operator()(bool output[], const S input[])
{
Vector<S, 2> _Point ;
_Point[0] = input[0];
_Point[1] = input[1]; //any point M
Vector<S, 2> _am=_Point – _PointA ;
Vector<S, 2> _bm=_Point – _PointB ;
Vector<S, 2> _cm=_Point – _PointC ;

S S_ABM = fabs(_ab[1]*_am[0]-_ab[0]*_am[1])/2. ; // the area of triangle ABM;
S S_BCM = fabs(_bc[1]*_bm[0]-_ca[0]*_bm[1])/2. ; // the area of triangle BCM;
S S_CAM = fabs(_ca[1]*_cm[0]-_ca[0]*_cm[1])/2. ; // the area of triangle CAM;
S Sum = fabs(_ca[1]*_ab[0]-_ca[0]*_ab[1])/2. ; // the area of triangle CAM;

output[0] = util::approxEqual(S_ABM + S_BCM + S_CAM, Sum) ;
return true;
}

template <typename S>
IndicatorTriangle2D<S>* createIndicatorTriangle2D(XMLreader const& params, bool verbose)
{
OstreamManager clout(std::cout,”createIndicatorTriangle2D”);
params.setWarningsOn(verbose);

Vector<S,2> PointA;
Vector<S,2> PointB;
Vector<S,2> PointC;

std::stringstream xmlCenter( params.getAttribute(“PointA”) );
xmlCenter >> PointA[0] >> PointA[1]>> PointB[0] >> PointB[1]>> PointC[0] >> PointC[1];

// std::stringstream xmlRadius( params.getAttribute(“length”) );
// xmlRadius >> xLength >> yLength;

return new IndicatorCuboid2D<S>(PointA, PointB, PointC);
}

[Adrian]

The next release is planned for sometime in the first half of 2022 (we realize that it already has been quite a while since the last release, however the next release will be worth the wait in more than one aspect… :-))

Your indicator class looks ok at first glance but I am afraid that we can not provide general C++ programming support here. Did you check it in a minimal example (i.e. without XML, as simple as you can make it). If you satisfy the basic indicator interface operator()(bool output[], const S input[]) then everything downstream should work.

[Fany]

That is great, really looking forward to the new release.

I change the way to creating the geometry, but I have a question of how to deal with BC. I found that the inclined boundary just had 3-4 voxels (but the resolution was actually 35), and the BC was determined by a thin inclined rectangle. I referred to the example “aorta3d”. as following:

[Author]

Posts