KNGeometricUnitCell.cpp

Go to the documentation of this file.

#include "CKNGlobal.h"
#include "KNGeometricUnitCell.h"
#include "KNGeometricAtomFactory.h"
#include "KNMPIManager.h"
#include "KNGeometricAtom.h"
#include <math.h>


double CKNGeometricUnitCell::m_fLength[3];
double CKNGeometricUnitCell::m_SubElementID = ATOM_DEFAULT_INDEX + 1;
int CKNGeometricUnitCell::m_nFullFillIndex = 15;

CKNGeometricUnitCell::CKNGeometricUnitCell()
{
    m_fUnitCellID = -1;
    m_bPeriodicElement = false;
    
    m_fAssignedIndex[_X] = 0;
    m_fAssignedIndex[_Y] = 0;
    m_fAssignedIndex[_Z] = 0;
}

CKNGeometricUnitCell::~CKNGeometricUnitCell()
{
    FinalSubElement();
}

void CKNGeometricUnitCell::SetLength(double fLength[3], double fOriginLength[3])
{ 
    m_fLength[_X] = fOriginLength[_X];
    m_fLength[_Y] = fOriginLength[_Y];
    m_fLength[_Z] = fOriginLength[_Z];
}

void CKNGeometricUnitCell::SetCoordination(CKNGeometricCoordination coordination)
{
    m_coordination.SetCoordination(coordination.GetCoordinationAll());
}

void CKNGeometricUnitCell::SetCoordination(double* pfCoordination)
{
    m_coordination.SetCoordination(pfCoordination);
}

void CKNGeometricUnitCell::SetCoordination(double fXCoordination, double fYCoordination, double fZCoordination)
{
    m_coordination.SetCoordination(fXCoordination, fYCoordination, fZCoordination);
}

void CKNGeometricUnitCell::ArrangeAtom(CKNCommandFileParser::LPINPUT_CMD_PARAM lpParam)
{
    int                             i, nSize;

    nSize = CKNGeometricAtomFactory::m_vectAtomList_1.size();
    for (i = 0; i < nSize; ++i)
    {
        CKNGeometricCoordination    AtomCoordnation = CKNGeometricAtomFactory::m_vectAtomList_1[i].GetCoordination();
        CKNGeometricAtom            SubElement;

        AtomCoordnation.SetOffset(m_coordination);
        SubElement.SetType(CKNGeometricAtom::A, CKNGeometricAtomFactory::GetMaterialNumber(lpParam->szDomainMat[0], CKNGeometricAtom::A));
        SubElement.SetDomainMaterial(CKNGeometricAtomFactory::GetMaterialNumber(lpParam->szDomainMat[0]));
        SubElement.SetCoordination(AtomCoordnation);
        m_vectIntenalAtom.push_back(SubElement);
    }

    nSize = CKNGeometricAtomFactory::m_vectAtomList_2.size();
    for (i = 0; i < nSize; ++i)
    {
        CKNGeometricCoordination    AtomCoordnation = CKNGeometricAtomFactory::m_vectAtomList_2[i].GetCoordination();
        CKNGeometricAtom            SubElement;

        AtomCoordnation.SetOffset(m_coordination);
        SubElement.SetCoordination(AtomCoordnation);
        SubElement.SetType(CKNGeometricAtom::C, CKNGeometricAtomFactory::GetMaterialNumber(lpParam->szDomainMat[0], CKNGeometricAtom::C));
        SubElement.SetDomainMaterial(CKNGeometricAtomFactory::GetMaterialNumber(lpParam->szDomainMat[0]));
        m_vectIntenalAtom.push_back(SubElement);
    }
}

void CKNGeometricUnitCell::SetSubDomainMaterial(CKNCommandFileParser::LPINPUT_CMD_PARAM lpParam)
{
    double              cx, cy, cz;
    int                 i, j, nSize = m_vectIntenalAtom.size();

    for (j = 0; j < lpParam->nSubDomainNumber; ++j)
    {
            for (i = 0; i < nSize; ++i)
            {
                CKNGeometricCoordination coordination = m_vectIntenalAtom[i].GetCoordination();
        
                cx = coordination.GetCoordination(_X);
                cy = coordination.GetCoordination(_Y);
                cz = coordination.GetCoordination(_Z);
                
                if (cx < lpParam->fOrigin[j+1][0] || cx > (lpParam->fOrigin[j+1][0] + lpParam->fShapeLength[j+1][0]))
                    continue;
        
                switch (lpParam->nShape[j + 1])
                {
                case BOX_SHAPE:
                    if (cy < lpParam->fOrigin[j + 1][1] || cy >(lpParam->fOrigin[j + 1][1] + lpParam->fShapeLength[j + 1][1]))
                        continue;

                    if (cz < lpParam->fOrigin[j + 1][2] || cz >(lpParam->fOrigin[j + 1][2] + lpParam->fShapeLength[j + 1][2]))
                        continue;
                    break;
                case CYLINDER_SHAPE:
                    if ((pow(((cy - lpParam->fOrigin[j + 1][_Y]) - (lpParam->fShapeLength[j + 1][_Y] / 2)), 2) / pow(lpParam->fShapeLength[j + 1][_Y] / 2, 2)) +
                        pow(((cz - lpParam->fOrigin[j + 1][_Z]) - (lpParam->fShapeLength[j + 1][_Z] / 2)), 2) / pow(lpParam->fShapeLength[j + 1][_Z] / 2, 2) > 1)
                        continue;
                    break;
                }

                MATERIAL_INDEX  nMaterialType;

                nMaterialType = CKNGeometricAtomFactory::GetMaterialNumber(lpParam->szDomainMat[j+1], m_vectIntenalAtom[i].GetType());
                m_vectIntenalAtom[i].SetMaterial(nMaterialType);
            }
    }

}

unsigned int CKNGeometricUnitCell::NumberingSubElement(CKNGeometricCoordination ShapeCoordination, double fLength[3], int nShapeType, char *pszDomainMatName)
{
    double              cx, cy, cz;
    int                 i, nSize = m_vectIntenalAtom.size();
    unsigned int        count = 0;

    for (i = 0; i < nSize; ++i)
    {
        CKNGeometricCoordination coordination = m_vectIntenalAtom[i].GetCoordination();
        
        cx = coordination.GetCoordination(_X);
        cy = coordination.GetCoordination(_Y);
        cz = coordination.GetCoordination(_Z);
                
        if (cx < 0 || cx > ShapeCoordination.GetCoordination(_X) + fLength[_X])
        {
            if( fabs(cx - (ShapeCoordination.GetCoordination(_X) + fLength[_X])) > GENERAL_TOLERANCE)
            {
                //printf("cx: %10.10f, cy: %10.10f, cz %10.10f\n", cx, cy, cz);
                continue;
            }
        }
        
        switch (nShapeType)
        {
        case BOX_SHAPE:
            if (cy < 0 || cy > ShapeCoordination.GetCoordination(_Y) + fLength[_Y])
                continue;

            if (cz < 0 || cz > ShapeCoordination.GetCoordination(_Z) + fLength[_Z])
                continue;
            break;
        case CYLINDER_SHAPE:
            if ((pow((cy - (fLength[_Y] / 2)), 2) / pow(fLength[_Y] / 2, 2)) +
                pow((cz - (fLength[_Z] / 2)), 2) / pow(fLength[_Z] / 2, 2) > 1)
                continue;
            break;
        }

        count++;

        MATERIAL_INDEX  nMaterialType;

        nMaterialType = CKNGeometricAtomFactory::GetMaterialNumber(pszDomainMatName, m_vectIntenalAtom[i].GetType());
        m_vectIntenalAtom[i].SetMaterial(nMaterialType);
        m_vectIntenalAtom[i].SetID(m_SubElementID++);
    }

    return count;
}

void CKNGeometricUnitCell::SetAssignIndex(double  fXAssignIndex, double fYAssignIndex, double fZAssignIndex)
{
    m_fAssignedIndex[_X] = floor(fXAssignIndex);
    m_fAssignedIndex[_Y] = floor(fYAssignIndex);
    m_fAssignedIndex[_Z] = floor(fZAssignIndex);
}

void CKNGeometricUnitCell::CheckingNeighborCandiate(double fAssignedCount[3], bool bFrontFace, bool bBackEnd)
{
    double              cx[3], cy[3], cz[3];
    int                 cxCount = 1, cyCount = 1, czCount = 1;
    int                 x, y, z;

    cx[0] = m_fAssignedIndex[_X];
    cy[0] = m_fAssignedIndex[_Y];
    cz[0] = m_fAssignedIndex[_Z];

    if (cx[0] > 0 || (cx[0] == 0 && false == bFrontFace))
        cx[cxCount++] = cx[0] - 1;

    if (cy[0] > 0)
        cy[cyCount++] = cy[0] - 1;

    if (cz[0] > 0)
        cz[czCount++] = cz[0] - 1;

    if (cx[0] < fAssignedCount[_X] - 1 || (cx[0] == fAssignedCount[_X] - 1 && false == bBackEnd))
        cx[cxCount++] = cx[0] + 1;

    if (cy[0] < fAssignedCount[_Y] - 1)
        cy[cyCount++] = cy[0] + 1;

    if (cz[0] < fAssignedCount[_Z] - 1)
        cz[czCount++] = cz[0] + 1;

    for (x = 0; x < cxCount; ++x)
    {
        for (y = 0; y < cyCount; ++y)
        {
            for (z = 0; z < czCount; ++z)
            {
                if (0 == x && 0 == y && 0 == z)
                    continue;

                m_vectNeighborCellID.push_back(CalculatingIndex(cx[x], cy[y], cz[z], fAssignedCount));
            }
        }
    }
}

CKNGeometricAtom* CKNGeometricUnitCell::GetSubElementAt(double fIndex)
{
    if (fIndex >= m_vectIntenalAtom.size())
        return NULL;

    return &m_vectIntenalAtom[(int)floor(fIndex)];
}

CKNGeometricAtom* CKNGeometricUnitCell::GetSubElementByID(double fID)
{
    int                         nSize, i;

    nSize = m_vectIntenalAtom.size();
    for( i = 0 ;i < nSize ; ++ i )
    {
        if( fID == m_vectIntenalAtom[i].GetID() )
            return &m_vectIntenalAtom[i];
    }

    return NULL;
}

void CKNGeometricUnitCell::FinalSubElement()
{
    m_vectIntenalAtom.clear();
}

double CKNGeometricUnitCell::CalculatingIndex(double fXAssignedIndex, double fYAssignedIndex, double fZAssignedIndex, double fAssignedCount[3])
{
    return fXAssignedIndex * floor(fAssignedCount[1]) * floor(fAssignedCount[2]) + fYAssignedIndex * floor(fAssignedCount[2]) + fZAssignedIndex;
}

void CKNGeometricUnitCell::CopyInnerAtomIndex(CKNGeometricUnitCell srcUnitCell)
{
    int             i, nSize = srcUnitCell.GetSubElementSize();
    for (i = 0; i < nSize; ++i)
        m_vectIntenalAtom[i].SetID(srcUnitCell.GetSubElementAt(i)->GetID());
}

void CKNGeometricUnitCell::ShiftAtomID(double fShift)
{
    int             i, nSize = m_vectIntenalAtom.size();

    for (i = 0; i < nSize; ++i)
        m_vectIntenalAtom[i].ShiftID(fShift);
}

void CKNGeometricUnitCell::ShiftID(double fShift)
{
    int             i, nSize = m_vectNeighborCellID.size();

    m_fUnitCellID += fShift; 
    for (i = 0; i < nSize; ++i)
        m_vectNeighborCellID[i] += fShift;
}

void CKNGeometricUnitCell::BuildNeighborInformation(double *pfAssignCount, std::vector<CKNGeometricAtom*> *pSurfaceList)
{
    double                                      i, j, nSecondIndex, nSize;
    int                                         nCheckFill = 0;
    int                                         nMatchingIndex;
    CKNGeometricAtom                            *pCompare1, *pCompare2, *pAtom;
    CKNGeometricAtom::ATOM_TYPE                 type1, type2;
    std::vector<CKNGeometricCoordination>       *pNeighborRelation = NULL;
    bool                                        bFindAll = false;

    nSize = m_vectIntenalAtom.size();
    for (i = 0; i < nSize; ++i)
    {
        if (ATOM_DEFAULT_INDEX == m_vectIntenalAtom[(unsigned int)i].GetID())
            continue;

        pCompare1 = &m_vectIntenalAtom[(unsigned int)i];
        type1 = pCompare1->GetType();

        if (CKNGeometricAtom::A == type1)
            pNeighborRelation = &CKNGeometricAtomFactory::m_vectNeighborRelation1To2;
        else
            pNeighborRelation = &CKNGeometricAtomFactory::m_vectNeighborRelation2To1;
        
        bFindAll = false;
        nCheckFill = 0;
        nSecondIndex = -1;
        for (j = 0; j < nSize - 1; ++j)
        {
            nSecondIndex++;

            if (nSecondIndex == i)
                nSecondIndex++;

            pCompare2 = &m_vectIntenalAtom[(unsigned int)nSecondIndex];

            type2 = pCompare2->GetType();
            if (type1 == type2)
                continue;

            nMatchingIndex = IsMachedNeighborRule(pCompare1, pCompare2, pNeighborRelation, nCheckFill);

            if (-1 != nMatchingIndex)
            {
                pCompare1->SetNeighbor(nMatchingIndex, pCompare2);
                nCheckFill |= (int)(pow((double)2, (int)nMatchingIndex));
            }

            if (m_nFullFillIndex == nCheckFill)
            {
                bFindAll = true;
                break;
            }
        }
                
        if (bFindAll)
            continue;

        int     nNeighborUnitcellSize = m_vectNeighborCellID.size();

        for (j = 0; j < nNeighborUnitcellSize; ++j)
        {
            double                      fNeighborID = m_vectNeighborCellID[(unsigned int)j];
            CKNGeometricUnitCell        NeighborCell;

#ifndef DISABLE_MPI_ROUTINE
            if (fNeighborID < 0)
            {
                if( m_vectUnitCellInPrev->size() > 0 )
                {
                    NeighborCell = m_vectUnitCellInPrev->at((int)fNeighborID + (int)(floor(pfAssignCount[_Y]) * floor(pfAssignCount[_Z])));
                    bFindAll = FindingNeighborInUnitcell(&NeighborCell, pCompare1, type1, pNeighborRelation, nCheckFill);
                }
            }
            else if (fNeighborID >= m_vectUnitCellInShape->size())
            {
                if( m_vectUnitCellInNext->size() > 0 )
                {
                    NeighborCell = m_vectUnitCellInNext->at((int)fNeighborID - (int)(floor(pfAssignCount[_X]) * floor(pfAssignCount[_Y]) * floor(pfAssignCount[_Z])));
                    bFindAll = FindingNeighborInUnitcell(&NeighborCell, pCompare1, type1, pNeighborRelation, nCheckFill);
                }
            }
            else
#endif //DISABLE_MPI_ROUTINE
            {
                NeighborCell = m_vectUnitCellInShape->at((unsigned int)fNeighborID);
                bFindAll = FindingNeighborInUnitcell(&NeighborCell, pCompare1, type1, pNeighborRelation, nCheckFill);
            }

            if (bFindAll)
                break;
        }
                
        if (bFindAll)
            continue;
    }

    for (i = 0; i < nSize; ++i)
    {
        int         nNeighborSize;
        if (ATOM_DEFAULT_INDEX == m_vectIntenalAtom[(unsigned int)i].GetID())
            continue;

        pAtom = &m_vectIntenalAtom[(unsigned int)i];
        nNeighborSize = pAtom->GetNeighborCount();
        for (j = 0; j < nNeighborSize; ++j)
        {
            if (ATOM_DEFAULT_INDEX == pAtom->GetNeighborAtomID(j))
            {
                pSurfaceList->push_back(pAtom);
                break;
            }
        }
    }
}

int CKNGeometricUnitCell::IsMachedNeighborRule(CKNGeometricAtom *pCompare1, CKNGeometricAtom *pCompare2, std::vector<CKNGeometricCoordination> *pNeighborRelation, int &nCheckFill)
{
    int                         i, nSize;
    CKNGeometricCoordination    coordination1;
    CKNGeometricCoordination    coordination2 = pCompare2->GetCoordination();
    
    nSize = pNeighborRelation->size();
    for (i = 0; i < nSize; ++i)
    {
        double      fNeighbor = pCompare1->GetNeighborAtomID(i);
        if (ATOM_DEFAULT_INDEX != fNeighbor)
            continue;

        coordination1 = pCompare1->GetCoordination();
        coordination1.SetOffset(pNeighborRelation->at(i));
        
        if (coordination2 == coordination1)
            return i;
    }

    return -1;
}

bool CKNGeometricUnitCell::FindingNeighborInUnitcell(CKNGeometricUnitCell *pUnitcell, CKNGeometricAtom *pCompare1, CKNGeometricAtom::ATOM_TYPE type1, std::vector<CKNGeometricCoordination> *pNeighborRelation, int &nCheckFill)
{
    int                                         i, nSize = pUnitcell->GetSubElementSize();
    CKNGeometricAtom::ATOM_TYPE                 type2;
    int                                         nMatchingIndex;

    for (i = 0; i < nSize; ++i)
    {
        CKNGeometricAtom *pCompare2 = pUnitcell->GetSubElementAt(i);
        type2 = pCompare2->GetType();
        /*if (type1 == type2)
            continue;*/

        nMatchingIndex = IsMachedNeighborRule(pCompare1, pCompare2, pNeighborRelation, nCheckFill);

        if (-1 != nMatchingIndex)
        {
            pCompare1->SetNeighbor(nMatchingIndex, pCompare2);
            nCheckFill |= (int)(pow((double)2, (int)nMatchingIndex));
        }

        if (15 == nCheckFill)
            return true;
    }

    return false;
}

void CKNGeometricUnitCell::Serialize(double *pBuffer)
{
    int             i, nSize = m_vectIntenalAtom.size();
    int             nPos = 0;

    for( i = 0; i < nSize;++i)
    {
        pBuffer[nPos++] = m_vectIntenalAtom[i].GetID();
        pBuffer[nPos++] = (double)m_vectIntenalAtom[i].GetType();
        pBuffer[nPos++] = (double)m_vectIntenalAtom[i].GetMaterialNumber();
        pBuffer[nPos++] = (double)m_vectIntenalAtom[i].IsPeriodicAtom();
        pBuffer[nPos++] = m_vectIntenalAtom[i].GetCoordination().GetCoordination(_X);
        pBuffer[nPos++] = m_vectIntenalAtom[i].GetCoordination().GetCoordination(_Y);
        pBuffer[nPos++] = m_vectIntenalAtom[i].GetCoordination().GetCoordination(_Z);
        pBuffer[nPos++] = (double)m_vectIntenalAtom[i].GetPeriodicDirection();
    }
}

void CKNGeometricUnitCell::Deserialize(double *pBuffer)
{
    int             i, nSize = CKNGeometricAtomFactory::GetAtomCountInUnitcell();
    int             nPos = 0;

    for( i = 0; i < nSize;++i)
    {
        CKNGeometricAtom        atom;

        atom.SetID(pBuffer[nPos++]);
        if( 0 == pBuffer[nPos++])
            atom.SetType(CKNGeometricAtom::A, (MATERIAL_INDEX)((int)pBuffer[nPos++]));
        else
            atom.SetType(CKNGeometricAtom::C, (MATERIAL_INDEX)((int)pBuffer[nPos++]));
        
        atom.SetPeriodic((bool)pBuffer[nPos++]);
        
        double      x = pBuffer[nPos++];
        double      y = pBuffer[nPos++];
        double      z = pBuffer[nPos++];
        atom.SetCoordination(x, y, z);
        atom.SetPeriodicDirection((int)pBuffer[nPos++]);
        m_vectIntenalAtom.push_back(atom);
    }
}

void CKNGeometricUnitCell::SetPeriodic(bool bPeriodic, int periodicDirection)
{
    int             i, nSize = CKNGeometricAtomFactory::GetAtomCountInUnitcell();
    int             nPos = 0;

    m_bPeriodicElement = bPeriodic;
    m_periodicDirection = periodicDirection;
    for (i = 0; i < nSize; ++i)
    {
        m_vectIntenalAtom[i].SetPeriodic(bPeriodic);
        m_vectIntenalAtom[i].SetPeriodicDirection(periodicDirection);
    }
}

void CKNGeometricUnitCell::SetUnitcellList(std::vector<CKNGeometricUnitCell> *pVect, std::vector<CKNGeometricUnitCell> *pVectPrev, std::vector<CKNGeometricUnitCell> *pVectNext)
{
    m_vectUnitCellInShape = pVect; 
    m_vectUnitCellInPrev = pVectPrev;
    m_vectUnitCellInNext = pVectNext;
}

void CKNGeometricUnitCell::SetPeriodicDirection(int periodicDirection)
{
    int             i, nSize = CKNGeometricAtomFactory::GetAtomCountInUnitcell();

    for (i = 0; i < nSize; ++i)
        m_vectIntenalAtom[i].SetPeriodicDirection(periodicDirection);
}

void CKNGeometricUnitCell::AddPeridoicDirection(CKNGeometricAtom::PERIODIC_DIRECTION periodicDirection)
{
    int             i, nSize = CKNGeometricAtomFactory::GetAtomCountInUnitcell();

    for (i = 0; i < nSize; ++i)
        m_vectIntenalAtom[i].AddPeridoicDirection(periodicDirection);
}