zAxis is upgraded and array2D is added for support for 2D array operations

- Cylinder geometric region is corrected for minPoint and maxPoint
This commit is contained in:
Hamidreza Norouzi 2024-04-17 14:04:34 -07:00
parent 79b987c711
commit c432880689
5 changed files with 253 additions and 158 deletions

View File

@ -63,6 +63,7 @@ containers/Field/Fields.cpp
containers/symArrayHD/symArrays.cpp
containers/List/anyList/anyList.cpp
structuredData/zAxis/zAxis.cpp
structuredData/box/box.cpp
structuredData/sphere/sphere.cpp
structuredData/cylinder/cylinder.cpp
@ -73,8 +74,6 @@ structuredData/peakableRegion/peakableRegion/peakableRegion.cpp
structuredData/pointStructure/internalPoints/internalPointsKernels.cpp
structuredData/pointStructure/internalPoints/internalPoints.cpp
structuredData/zAxis/zAxis.cpp
structuredData/line/line.cpp
structuredData/infinitePlane/infinitePlane.cpp
structuredData/plane/plane.cpp

View File

@ -27,7 +27,6 @@ FUNCTION_H
bool pFlow::cylinder::calculateParams()
{
WARNING<<"Use of cylinder requires modifications to zAxis"<<END_WARNING;
auto p1p2 = p2_ - p1_;
if( p1p2.length() > smallValue )
@ -45,8 +44,11 @@ bool pFlow::cylinder::calculateParams()
realx3 minPinZ(-sqrt(radius2_), -sqrt(radius2_), 0.0);
realx3 maxPinZ( sqrt(radius2_), sqrt(radius2_), sqrt(axisVector2_));
minPoint_ = zA.transferBackZ(minPinZ);
maxPoint_ = zA.transferBackZ(maxPinZ);
auto minp = zA.transferBackZ(minPinZ);
auto maxp = zA.transferBackZ(maxPinZ);
minPoint_ = min(minp, maxp);
maxPoint_ = max(minp, maxp);
return true;
}

View File

@ -0,0 +1,188 @@
#ifndef __array2D_hpp__
#define __array2D_hpp__
#include "iOstream.hpp"
namespace pFlow
{
template<typename T, size_t nRow, size_t nCol>
struct array2D
{
T elements_[nRow][nCol];
inline
T& operator()(size_t i, size_t j) noexcept
{
return elements_[i][j];
}
inline
const T& operator()(size_t i, size_t j)const noexcept
{
return elements_[i][j];
}
constexpr size_t size()const noexcept
{
return nRow*nCol;
}
constexpr size_t nCols()const noexcept
{
return nCols;
}
constexpr size_t nRows()const noexcept
{
return nRows;
}
};
template<typename T, size_t nRow, size_t nCol>
array2D<T, nRow, nCol> operator+(
const array2D<T, nRow, nCol>& arr1,
const array2D<T, nRow, nCol>& arr2)
{
array2D<T, nRow, nCol> res;
for(size_t i=0; i<nRow; i++)
for(size_t j=0; j<nCol; j++)
res(i,j)= arr1(i,j)+arr2(i,j);
return res;
}
template<typename T, size_t nRow, size_t nCol>
array2D<T, nRow, nCol> operator-(
const array2D<T, nRow, nCol>& arr1,
const array2D<T, nRow, nCol>& arr2)
{
array2D<T, nRow, nCol> res;
for(size_t i=0; i<nRow; i++)
for(size_t j=0; j<nCol; j++)
res(i,j)= arr1(i,j)-arr2(i,j);
return res;
}
template<typename T, size_t nRow, size_t nCol>
array2D<T, nRow, nCol> operator*(
const array2D<T, nRow, nCol>& arr1,
const array2D<T, nRow, nCol>& arr2)
{
array2D<T, nRow, nCol> res;
for(size_t i=0; i<nRow; i++)
for(size_t j=0; j<nCol; j++)
res(i,j)= arr1(i,j)*arr2(i,j);
return res;
}
template<typename T, size_t nRow, size_t nCol>
array2D<T, nRow, nCol> operator*(
const T& s,
const array2D<T, nRow, nCol>& arr2)
{
array2D<T, nRow, nCol> res;
for(size_t i=0; i<nRow; i++)
for(size_t j=0; j<nCol; j++)
res(i,j)= s*arr2(i,j);
return res;
}
template<typename T, size_t nRow, size_t nCol>
array2D<T, nRow, nCol> operator/(
const array2D<T, nRow, nCol>& arr1,
const array2D<T, nRow, nCol>& arr2)
{
array2D<T, nRow, nCol> res;
for(size_t i=0; i<nRow; i++)
for(size_t j=0; j<nCol; j++)
res(i,j)= arr1(i,j)/arr2(i,j);
return res;
}
template<typename T, size_t nRow, size_t nCol>
array2D<T, nRow, nCol> operator/(
const T& val,
const array2D<T, nRow, nCol>& arr2)
{
array2D<T, nRow, nCol> res;
for(size_t i=0; i<nRow; i++)
for(size_t j=0; j<nCol; j++)
res(i,j)= val/arr2(i,j);
return res;
}
template <typename T, size_t nRow, size_t nInner, size_t nCol >
array2D<T, nRow, nCol>
MatMul(
const array2D<T, nRow, nInner>& A,
const array2D<T, nInner, nCol>& B)
{
array2D<T, nRow, nCol> C;
for (size_t row = 0; row < nRow; row++)
{
for (size_t col = 0; col < nCol; col++)
{
T sum = 0;
for (size_t inner = 0; inner < nInner; inner++)
{
sum += A(row,inner) * B(inner,col);
}
C(row,col) = sum;
}
}
return C;
}
template<typename T, size_t nRow, size_t nCol>
array2D<T, nCol, nRow>
transpose(const array2D<T, nRow, nCol>& arr)
{
array2D<T, nCol, nRow> tArr;
for(size_t i=0; i<nRow; i++)
{
for(size_t j=0; j<nCol; j++)
{
tArr(j,i) = arr(i,j);
}
}
return tArr;
}
template<typename T, size_t nRow, size_t nCol>
iOstream& operator<<(iOstream& os, const array2D<T, nRow, nCol>& arr)
{
os<<'[';
for(size_t i=0; i<nRow; i++)
{
os<<'[';
for(size_t j=0; j<nCol-1; j++)
{
os<<arr(i,j)<<' ';
}
if(i < nRow-1)
os<<arr(i,nCol-1)<<"]\n";
else
os<<arr(i,nCol-1)<<"]";
}
os<<']';
return os;
}
}
#endif //__array2D_hpp__

View File

@ -22,129 +22,60 @@ Licence:
#include "zAxis.hpp"
pFlow::zAxis::zAxis(const realx3 &p1, const realx3 &p2) :
p1_(p1),
p2_(p2)
pFlow::zAxis::zAxis(const realx3 &p1, const realx3 &p2)
:
p1_(p1),
p2_(p2)
{
n_ = p2-p1;
auto len = pFlow::length(n_);
if(len < smallValue )
if( equal(p1,p2))
{
fatalErrorInFunction<<
"points are equal "<< p1 <<" "<<p2<<endl;
fatalExit;
}
n_ /= len;
makeTransMatrix();
makeTransMatrix(p2-p1);
}
pFlow::realx3 pFlow::zAxis::transferToZ(const realx3 & p)
pFlow::realx3 pFlow::zAxis::transferToZ(const realx3 & p)const
{
real pp[4][1] = {p.x(), p.y(), p.z(), 1.0};
real pn[4][1];
const auto p0 = p1_-p;
array2D<real,3,1> po {p0.x(), p0.y(), p0.z()};
MatMul(Trans_z_xz_P1_, pp, pn);
auto tp = MatMul(rotMat_, po);
return realx3(
pn[0][0],
pn[1][0],
pn[2][0]
);
return realx3(tp(0,0), tp(1,0), tp(2,0));
}
pFlow::realx3 pFlow::zAxis::transferBackZ(const realx3 & p)
{
real pp[4][1] = { p.x(), p.y(), p.z(), 1.0 };
real pn[4][1];
MatMul(ITrans_P1_xz_z_, pp, pn);
return realx3(
pn[0][0],
pn[1][0],
pn[2][0]
);
}
void pFlow::zAxis::makeTransMatrix()
pFlow::realx3 pFlow::zAxis::transferBackZ(const realx3 & p)const
{
// transfering point p1 to the origin
real TransP1[4][4] =
{
1.0, 0.0, 0.0, -p1_.x(),
0.0, 1.0, 0.0, -p1_.y(),
0.0, 0.0, 1.0, -p1_.z(),
0.0, 0.0, 0.0, 1.0
};
auto rp = MatMul(invRotMat_,array2D<real,3,1>{p.x(), p.y(), p.z()});
// for transformation back
real ITransP1[4][4] =
{
1.0, 0.0, 0.0, p1_.x(),
0.0, 1.0, 0.0, p1_.y(),
0.0, 0.0, 1.0, p1_.z(),
0.0, 0.0, 0.0, 1.0
};
return realx3(rp(0,0)+p1_.x(),rp(1,0)+p1_.y(), rp(2,0)+p1_.z());
}
real u = n_.x();
real v = n_.y();
real w = n_.z();
real u2v2 = sqrt(u*u + v*v);
//correcting the transformation matrix in the case of coincidence with z - axis
if ( equal(w,1.0) )
{
assignMat(TransP1 , Trans_z_xz_P1_);
assignMat(ITransP1, ITrans_P1_xz_z_);
return;
}
u2v2 = max(smallValue, u2v2);
real TransXZ[4][4] =
{
u / u2v2, v / u2v2, 0.0, 0.0,
-v / u2v2, u / u2v2, 0.0, 0.0,
0.0, 0.0, 1.0, 0.0,
0.0, 0.0, 0.0, 1.0
};
real TransZ[4][4] =
{
w, 0.0, -u2v2, 0.0,
0.0, 1.0, 0.0, 0.0,
u2v2, 0.0, w, 0.0,
0.0, 0.0, 0.0, 1.0
};
real temp[4][4];
// creat transformation matrix to transfer point from line axis to z-axis
MatMul(TransXZ, TransP1, temp);
MatMul(TransZ, temp, Trans_z_xz_P1_);
real ITransXZ[4][4] =
{
u / u2v2, -v / u2v2, 0.0, 0.0,
+v / u2v2, u / u2v2, 0.0, 0.0,
0.0, 0.0, 1.0, 0.0,
0.0, 0.0, 0.0, 1.0
};
void pFlow::zAxis::makeTransMatrix(const realx3& vec)
{
const realx3 Unit = {0,0,1};
const auto v_unit = normalize(vec);
auto uvw = cross(v_unit,Unit);
real ITransZ[4][4] =
{
w, 0.0, +u2v2, 0.0,
0.0, 1.0, 0.0, 0.0,
-u2v2, 0.0, w, 0.0,
0.0, 0.0, 0.0, 1.0
};
const auto rcos = dot(v_unit, Unit);
const auto rsin = uvw.length();
// creat transformation matrix to transfer point to from z-axis to line axis
MatMul(ITransXZ, ITransZ, temp);
MatMul(ITransP1, temp, ITrans_P1_xz_z_);
if (abs(rsin) > smallValue)
uvw /= rsin;
const auto& [u, v, w] = uvw;
rotMat_ =
ArrayType({rcos, 0, 0, 0, rcos, 0, 0, 0, rcos})+
ArrayType({0, -w, v, w, 0, -u, -v, u, 0})+
(1-rcos)*ArrayType({u*u, u*v, u*w, u*v, v*v, w*v, u*w, v*w, w*w});
invRotMat_ = transpose(rotMat_);
}

View File

@ -29,78 +29,53 @@ p1 and p2
#define __zAxis_hpp__
#include "types.hpp"
#include "array2D.hpp"
namespace pFlow
{
template <typename T, int32 nRow, int32 nInner, int32 nCol >
void MatMul(T(&A)[nRow][nInner], T(&B)[nInner][nCol], T(&C)[nRow][nCol]);
template <typename T, int32 nRow, int32 nCol >
void assignMat(T(&A)[nRow][nCol], T(&B)[nRow][nCol]);
class zAxis
{
private:
using ArrayType = array2D<real,3uL,3uL>;
/// the origin of the transformed coord
realx3 p1_;
/// the direction vector of rotated coordinates
realx3 p2_;
/// transformation matrix to rotate original coordinates
/// to rotated coordinates
ArrayType rotMat_;
/// rotation matrix to rotate back from rotated coordinates
/// to the original axis
ArrayType invRotMat_;
void makeTransMatrix(const realx3& v);
public:
// constructors
zAxis(const realx3 &lp1, const realx3 &lp2);
inline
real length()const
{
return pFlow::length(p2_-p1_);
}
realx3 transferToZ(const realx3 & p);
realx3 transferToZ(const realx3 & p)const;
realx3 transferBackZ(const realx3 & p);
realx3 transferBackZ(const realx3 & p)const;
private:
void makeTransMatrix();
protected:
realx3 p1_;
realx3 p2_;
realx3 n_;
real Trans_z_xz_P1_[4][4];
real ITrans_P1_xz_z_[4][4];
};
template <typename T, int32 nRow, int32 nInner, int32 nCol >
void MatMul(T(&A)[nRow][nInner], T(&B)[nInner][nCol], T(&C)[nRow][nCol])
{
for (int32 row = 0; row < nRow; row++)
{
for (int32 col = 0; col < nCol; col++)
{
T sum = 0;
for (int inner = 0; inner < nInner; inner++)
{
sum += A[row][inner] * B[inner][col];
}
C[row][col] = sum;
}
}
}
template <typename T, int32 nRow, int32 nCol >
void assignMat(T(&A)[nRow][nCol], T(&B)[nRow][nCol])
{
for (int32 row = 0; row < nRow; row++)
{
for (int32 col = 0; col < nCol; col++)
{
B[row][col] = A[row][col];
}
}
}
}
#endif