RCS Computation Toolkit This repository provides C++ and Python tools for computing radar cross section (RCS) from surface currents or analytical models
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 

640 lines
16 KiB

#ifndef GROUP_H
#define GROUP_H
#include "vtr.h"
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <set>
using namespace std;
template<class T>
class Basis {
private:
int ID;
int nTri;
T **triPtr;
int cnt;
public:
// set procedures
Basis(){
ID = -1;
nTri = 0;
triPtr = NULL;
cnt = 0;
}
void setID(int id) { ID = id; }
void InitialtriPtr(int n){
if (triPtr != NULL)
delete [] triPtr;
triPtr = new T* [n];
}
void AddTriPtr(T *tt) {
int i;
for (i = 0; i < cnt; i ++){
if (triPtr[i] == tt)
return;
}
triPtr[cnt] = tt;
cnt ++;
}
void AddnTri() { nTri ++; }
// Get procedures
int getID() { return ID; }
int getnTri() { return nTri; }
int getCnt() { return cnt; }
T* getTriPtr(int i) { return triPtr[i]; }
};
typedef struct box{
float xmin, ymin, zmin;
float xmax, ymax, zmax;
vtr center;
} box;
// element: tri3d or face or tetra
// In this subroutine, I used data struture of tri3d
template<class element>
class group {
private:
box BoundBox;
group<element> *children[2];
group<element> *mother;
element **elemArray;
int nElement, nBasis, nChild;
int *bArray, counter;
int IdMax, IdMin;
double radius;
bool basisBuilt;
bool PointInBox(double, double, double);
void EqualPartition(int, double, double, double &);
void AddnElement();
public:
group();
void MakeChildren();
void FormBox();
void FormBoundary(double, double, double, double, double, double);
void Initialize();
void InsertElement(element *);
void ComputeGroupRadius();
void PartitionMultiLevelSVD(int, int&);
void BuildGroupBarray();
void BuildMotherBarray();
void findGroup(group<element> **);
void findGroupHelper(int &, group<element> **, group<element> *);
// set functions
void SetnBasis(int n ) { nBasis = n; }
void SetbArray(int *array) { bArray = array; }
void SetnElement(int n) { nElement = n; }
void SetMother(group<element> *grp) { mother = grp; }
void SetMaxId(int n) { IdMax = n; }
void SetMinId(int n) { IdMin = n; }
void SetBasisBuilt() { basisBuilt = true; }
//get functions
box GetBox() { return BoundBox; }
group *GetMother() { return mother; }
group **GetChildren() { return children; }
group *GetChildren(int i) { return children[i]; }
vtr GetGroupCenter() { return BoundBox.center; }
element* GetElemArray(int i) { return elemArray[i]; }
int GetnElement() { return nElement; }
int GetnChild() { return nChild; }
int GetnBasis() { return nBasis; }
int* GetbArray() { return bArray; }
double GetRadius() { return radius; }
int GetMaxId() { return IdMax; }
int GetMinId() { return IdMin; }
bool IsBasisBuilt() { return basisBuilt; }
};
template<class element>
group<element>::group()
{
basisBuilt = false;
nElement = 0;
counter = 0;
nBasis = 0;
IdMin = int (1.0e6);
IdMax = -1;
elemArray = 0;
nChild = 2;
mother = 0;
for (int i = 0; i < nChild; i ++)
children[i] = 0;
}
template<class element>
void group<element>::FormBoundary(double x_min, double x_max,
double y_min, double y_max,
double z_min, double z_max)
{
BoundBox.xmin = x_min; BoundBox.xmax = x_max;
BoundBox.ymin = y_min; BoundBox.ymax = y_max;
BoundBox.zmin = z_min; BoundBox.zmax = z_max;
}
template<class element>
void group<element>::FormBox()
{
double x, y, z;
BoundBox.xmin = 1.0e10;
BoundBox.ymin = 1.0e10;
BoundBox.zmin = 1.0e10;
BoundBox.xmax = -1.0e10;
BoundBox.ymax = -1.0e10;
BoundBox.zmax = -1.0e10;
int i, j;
for (i = 0; i < nElement; i ++){
element *elem = elemArray[i];
for (j = 0 ; j < 3; j ++){
x = elem->GetNodePtr(j)->GetX();
y = elem->GetNodePtr(j)->GetY();
z = elem->GetNodePtr(j)->GetZ();
if (x < BoundBox.xmin) BoundBox.xmin = x;
if (y < BoundBox.ymin) BoundBox.ymin = y;
if (z < BoundBox.zmin) BoundBox.zmin = z;
if (x > BoundBox.xmax) BoundBox.xmax = x;
if (y > BoundBox.ymax) BoundBox.ymax = y;
if (z > BoundBox.zmax) BoundBox.zmax = z;
}
}
x = (BoundBox.xmin + BoundBox.xmax) / 2.0;
y = (BoundBox.ymin + BoundBox.ymax) / 2.0;
z = (BoundBox.zmin + BoundBox.zmax) / 2.0;
BoundBox.center.setvtr(x, y, z);
// compute group radius
vtr Gc = BoundBox.center;
vtr point[8];
radius = 0.0;
double xmax, xmin, ymax, ymin, zmax, zmin;
xmax = BoundBox.xmax;
ymax = BoundBox.ymax;
zmax = BoundBox.zmax;
xmin = BoundBox.xmin;
ymin = BoundBox.ymin;
zmin = BoundBox.zmin;
point[0].setvtr(xmin, ymin, zmin);
point[1].setvtr(xmin, ymax, zmin);
point[2].setvtr(xmax, ymin, zmin);
point[3].setvtr(xmax, ymax, zmin);
point[4].setvtr(xmin, ymin, zmax);
point[5].setvtr(xmin, ymax, zmax);
point[6].setvtr(xmax, ymin, zmax);
point[7].setvtr(xmax, ymax, zmax);
for (i = 0 ; i < 8 ; i ++) {
double length = Length(point[i] - Gc);
if (length > radius)
radius = length;
}
}
/* // update the coordinate for the bounding box, */
/* // as well as the center of the box */
/* template<class element> */
/* void group<element>::FormBox() */
/* { */
/* double x, y, z; */
/* BoundBox.xmin = 1.0e10; */
/* BoundBox.ymin = 1.0e10; */
/* BoundBox.zmin = 1.0e10; */
/* BoundBox.xmax = -1.0e10; */
/* BoundBox.ymax = -1.0e10; */
/* BoundBox.zmax = -1.0e10; */
/* int i; */
/* for (i = 0; i < nElement; i ++){ */
/* element *elem = elemArray[i]; */
/* elem->GetCenter(x, y, z); */
/* if (x < BoundBox.xmin) BoundBox.xmin = x; */
/* if (y < BoundBox.ymin) BoundBox.ymin = y; */
/* if (z < BoundBox.zmin) BoundBox.zmin = z; */
/* if (x > BoundBox.xmax) BoundBox.xmax = x; */
/* if (y > BoundBox.ymax) BoundBox.ymax = y; */
/* if (z > BoundBox.zmax) BoundBox.zmax = z; */
/* } */
/* x = (BoundBox.xmin + BoundBox.xmax) / 2.0; */
/* y = (BoundBox.ymin + BoundBox.ymax) / 2.0; */
/* z = (BoundBox.zmin + BoundBox.zmax) / 2.0; */
/* BoundBox.center.setvtr(x, y, z); */
/* // compute group radius */
/* vtr Gc = BoundBox.center; */
/* vtr point[8]; */
/* radius = 0.0; */
/* double xmax, xmin, ymax, ymin, zmax, zmin; */
/* xmax = BoundBox.xmax; */
/* ymax = BoundBox.ymax; */
/* zmax = BoundBox.zmax; */
/* xmin = BoundBox.xmin; */
/* ymin = BoundBox.ymin; */
/* zmin = BoundBox.zmin; */
/* point[0].setvtr(xmin, ymin, zmin); */
/* point[1].setvtr(xmin, ymax, zmin); */
/* point[2].setvtr(xmax, ymin, zmin); */
/* point[3].setvtr(xmax, ymax, zmin); */
/* point[4].setvtr(xmin, ymin, zmax); */
/* point[5].setvtr(xmin, ymax, zmax); */
/* point[6].setvtr(xmax, ymin, zmax); */
/* point[7].setvtr(xmax, ymax, zmax); */
/* for (i = 0 ; i < 8 ; i ++) { */
/* double length = Length(point[i] - Gc); */
/* if (length > radius) */
/* radius = length; */
/* } */
/* } */
// initialize the elemArray
template<class element>
void group<element>::Initialize()
{
if (elemArray != 0) delete [] elemArray;
elemArray = new element* [nElement];
}
// increment the number of element
template<class element>
void group<element>::AddnElement()
{
nElement ++;
}
// insert an element to elemArray
template<class element>
void group<element>::InsertElement(element *T)
{
/*
int i;
for (i = 0; i < nElement; i ++) {
if (elemArray[i] == T)
return;
}
*/
elemArray[counter ++] = T;
}
// Compute the radius of bounding box
template<class element>
void group<element>::ComputeGroupRadius()
{
vtr Gc = BoundBox.center;
vtr point[8];
radius = 0.0;
double xmax, xmin, ymax, ymin, zmax, zmin;
xmax = BoundBox.xmax;
ymax = BoundBox.ymax;
zmax = BoundBox.zmax;
xmin = BoundBox.xmin;
ymin = BoundBox.ymin;
zmin = BoundBox.zmin;
point[0].setvtr(xmin, ymin, zmin);
point[1].setvtr(xmin, ymax, zmin);
point[2].setvtr(xmax, ymin, zmin);
point[3].setvtr(xmax, ymax, zmin);
point[4].setvtr(xmin, ymin, zmax);
point[5].setvtr(xmin, ymax, zmax);
point[6].setvtr(xmax, ymin, zmax);
point[7].setvtr(xmax, ymax, zmax);
int i;
for (i = 0 ; i < 8 ; i ++) {
double length = Length(point[i] - Gc);
if (length > radius)
radius = length;
}
}
// determine if a point inside the bounding box
template<class element>
bool group<element>::PointInBox(double x, double y, double z)
{
double xmax, xmin, ymax, ymin, zmax, zmin;
xmax = BoundBox.xmax;
ymax = BoundBox.ymax;
zmax = BoundBox.zmax;
xmin = BoundBox.xmin;
ymin = BoundBox.ymin;
zmin = BoundBox.zmin;
bool In = true, Out = false;
if (xmin <= x && xmax >= x &&
ymin <= y && ymax >= y
&& zmin <= z && zmax >= z)
return In;
return Out;
}
// this subroutine determines where to cut the box into half
template<class element>
void group<element>::EqualPartition(int index, double min,
double max, double &mid)
{
int i;
double center, xc, yc, zc;
int check = 0;
int leftcnt;
bool found = true;
float range_min = 0.45;
float range_max = 0.55;
double midMax = max;
double midMin = min;
while (found) {
if (check == 1) break;
leftcnt = 0;
for (i = 0; i < nElement; i ++) {
element *elem = elemArray[i];
elem->GetCenter(xc, yc, zc);
if (index == 0) center = xc;
if (index == 1) center = yc;
if (index == 2) center = zc;
if (check == 0) mid = 0.5 * (max + min);
if ((min <= center) && (mid >= center)) {
leftcnt ++;
}
}
double Cnt = leftcnt;
if ((Cnt >= range_min * nElement) && (Cnt <= range_max * nElement))
break;
if (Cnt <= range_min * nElement) {
midMin = mid;
mid = 0.5 * (midMax + mid);
}
if (Cnt >= range_max * nElement) {
midMax = mid;
mid = 0.5 * (mid + midMin);
}
check ++;
}
}
template<class element>
void group<element>::MakeChildren()
{
int i, j;
for (i = 0; i < nChild; i ++)
children[i] = new group<element> [1];
double xmax, xmin, ymax, ymin, zmax, zmin;
xmax = BoundBox.xmax;
ymax = BoundBox.ymax;
zmax = BoundBox.zmax;
xmin = BoundBox.xmin;
ymin = BoundBox.ymin;
zmin = BoundBox.zmin;
double xmid, ymid, zmid;
double dmax = xmax - xmin;
int optCASE = 0;
if ((ymax - ymin) > dmax) {
optCASE = 1;
dmax = ymax - ymin;
}
if ((zmax - zmin) > dmax) {
optCASE = 2;
dmax = zmax - zmin;
}
if (dmax < 1.0e-8) return;
if (optCASE == 0) {
EqualPartition(0, xmin, xmax, xmid);
children[0]->FormBoundary(xmin, xmid, ymin, ymax, zmin, zmax);
children[1]->FormBoundary(xmid, xmax, ymin, ymax, zmin, zmax);
} else if (optCASE == 1) {
EqualPartition(1, ymin, ymax, ymid);
children[0]->FormBoundary(xmin, xmax, ymin, ymid, zmin, zmax);
children[1]->FormBoundary(xmin, xmax, ymid, ymax, zmin, zmax);
} else {
EqualPartition(2, zmin, zmax, zmid);
children[0]->FormBoundary(xmin, xmax, ymin, ymax, zmin, zmid);
children[1]->FormBoundary(xmin, xmax, ymin, ymax, zmid, zmax);
}
double xc, yc, zc;
for (i = 0; i < nElement; i ++) {
element *elem = elemArray[i];
elem->GetCenter(xc, yc, zc);
if (fabs(xmin - xmax) < 1.0e-8) xc = xmin;
if (fabs(ymin - ymax) < 1.0e-8) yc = ymin;
if (fabs(zmin - zmax) < 1.0e-8) zc = zmin;
bool InBox;
InBox = children[0]->PointInBox(xc, yc, zc);
if (InBox) {
children[0]->AddnElement();
} else {
children[1]->AddnElement();
}
}
int nElem = 0;
for (i = 0; i < nChild; i ++) {
int elem = children[i]->GetnElement();
if (elem == 0) {
/* cout << "one of children is empty\n"; */
for (j = 0; j < nChild; j ++)
children[j] = 0;
return;
}
nElem += elem;
}
if (nElem != nElement) {
cout << "Error : Mother is not splitted into children" << endl;
cout << nElem << " " << nElement << endl;
exit(1);
}
for (i = 0; i < nChild; i ++) {
children[i]->SetMother(this);
children[i]->Initialize();
}
set<element*> child0;
child0.clear();
set<element*> child1;
child1.clear();
for (i = 0; i < nElement; i ++) {
element *elem = elemArray[i];
elem->GetCenter(xc, yc, zc);
if (fabs(xmin - xmax) < 1.0e-8) xc = xmin;
if (fabs(ymin - ymax) < 1.0e-8) yc = ymin;
if (fabs(zmin - zmax) < 1.0e-8) zc = zmin;
bool InBox;
InBox = children[0]->PointInBox(xc, yc, zc);
if (InBox) {
// children[0]->InsertElement(elem);
child0.insert(elem);
} else {
// children[1]->InsertElement(elem);
child1.insert(elem);
}
}
int cntt=0;
for( typename set<element*>::iterator it= child0.begin() ; it!=child0.end(); it++){
children[0]->elemArray[cntt]= *it;
cntt++;
}
cntt=0;
for(typename set<element*>::iterator it=child1.begin();it!=child1.end();it++){
children[1]->elemArray[cntt]= *it;
cntt++;
}
child0.clear();
child1.clear();
}
template<class element>
void group<element>::PartitionMultiLevelSVD(int grpsize, int &numGrp)
{
int i;
if (nElement <= grpsize){
if (nElement != 0) {
FormBox();
numGrp ++;
}
return;
}
FormBox();
MakeChildren();
int noChild = 0;
for (i = 0; i < nChild; i ++) {
if (children[i] != 0)
children[i]->PartitionMultiLevelSVD(grpsize, numGrp);
else noChild = 1;
}
if (noChild && nElement != 0) {
FormBox();
numGrp ++;
return;
}
}
template<class element>
void group<element>::BuildGroupBarray()
{
int i, j;
for (i = 0; i < nElement; i ++) {
element *elem = elemArray[i];
int edgeNum = 3;
for (j = 0; j < edgeNum; j ++) {
int id = elem->GetEdgePtr(j)->GetId();
int VALID = elem->GetEdgePtr(j)->GetNegativeId();
if (id < 0 || VALID != -1) continue;
elem->GetEdgePtr(j)->SetNegativeId(-2);
nBasis ++;
if (id < IdMin) IdMin = id;
if (id > IdMax) IdMax = id;
}
}
bArray = new int [nBasis];
int cnt = 0;
for (i = 0; i < nElement; i ++) {
element *elem = elemArray[i];
int edgeNum = 3;
for (j = 0; j < edgeNum; j ++) {
int id = elem->GetEdgePtr(j)->GetId();
int VALID = elem->GetEdgePtr(j)->GetNegativeId();
if (id < 0 || VALID != -2) continue;
elem->GetEdgePtr(j)->SetNegativeId(1);
bArray[id - IdMin] = id;
cnt ++;
}
}
basisBuilt = true;
}
// Combine the children unknown IDs into mother
template<class element>
void group<element>::BuildMotherBarray()
{
if (mother == 0) return;
bool built = mother->IsBasisBuilt();
if ( built ) return;
group<element> **grp = mother->GetChildren();
int n1 = grp[0]->GetnBasis();
int n2 = grp[1]->GetnBasis();
if (n1 == 0 || n2 == 0) return;
int *ID1 = grp[0]->GetbArray();
int *ID2 = grp[1]->GetbArray();
int n = n1 + n2;
mother->SetnBasis(n);
int *array = new int [n];
int min1 = grp[0]->GetMinId();
int min2 = grp[1]->GetMinId();
int min = (min1 < min2) ? min1 : min2;
int max1 = grp[0]->GetMaxId();
int max2 = grp[1]->GetMaxId();
int max = (max1 > max2) ? max1 : max2;
int i;
for (i = 0; i < n1; i ++)
array[ID1[i] - min] = ID1[i];
for (i = 0; i < n2; i ++)
array[ID2[i] - min] = ID2[i];
mother->SetbArray(array);
mother->SetBasisBuilt();
mother->SetMinId(min);
mother->SetMaxId(max);
mother->BuildMotherBarray();
}
template<class element>
void group<element>::findGroup(group<element> **grp)
{
int cnt = 0;
findGroupHelper(cnt, grp, this);
}
template<class element>
void group<element>::findGroupHelper(int &cnt, group<element> **grp,
group<element> *mgp)
{
int flag = 0;
if (mgp->children[0] == 0) flag = 1;
else
findGroupHelper(cnt, grp, mgp->children[0]);
if (mgp->children[1] == 0) {
if (flag == 1 && mgp->GetnElement() != 0)
grp[cnt++] = mgp;
}
else findGroupHelper(cnt, grp, mgp->children[1]);
}
#endif