/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Application
    foamMeshToAbaqus

Description
    Converts an OpenFOAM mesh to an Abaqus input file.
    Creates a node set and and element set and a surface for each boundary
    patch.
    Also creates a element set for each material in the materials file (if
    it is exists).
    Only works for hexahedral cells as yet.
    Checked for Abaqus-6.9-2.

Author
    philip.cardiff@ucd.ie

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "OFstream.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
#   include "setRootCase.H"
#   include "createTime.H"
#   include "createMesh.H"

    //----------------------------------------------------------------------//
    //- open abaqus input file
    //----------------------------------------------------------------------//
    Info << "Opening Abaqus input file" << endl << endl;
    OFstream abaqusInp("abaqusMesh.inp");

    //- input header
    Info << "Writing input file header" << endl << endl;
    abaqusInp << "*Heading\n"
        << "** Job name: OpenFoamMesh Model name: OpenFoamMesh\n"
        << "** Generated by: Abaqus/CAE 6.9-2\n"
        << "*Preprint, echo=NO, model=NO, history=NO, contact=NO\n"
        << "**\n"
        << "** PARTS\n"
        << "**\n"
        << "*Part, name=OpenFoamMeshPart"
        << endl;


    //----------------------------------------------------------------------//
    //- write nodes
    //----------------------------------------------------------------------//
    Info << "Writing Nodes" << endl << endl;
    abaqusInp << "*Node" << endl;
    const pointField& points = mesh.points();
    forAll(points, pointi)
    {
        abaqusInp << "\t" << (pointi+1)
            << ",\t" << (points[pointi]).x()
            << ",\t" << (points[pointi]).y()
            << ",\t" << (points[pointi]).z()
            << endl;
    }


    //----------------------------------------------------------------------//
    //- determine abaqusCellPoints
    //----------------------------------------------------------------------//
    //- for hex 1st order elements, abaqus orders the points
    //- where nodes 1 2 3 4 are clockwise from the outside of the cell
    //- and then 5 6 7 8 are apposite 1 2 3 4 respectively
    Info << "Determining Abaqus element node ordering" << endl << endl;

    const cellList& cells = mesh.cells();
    const faceList& faces = mesh.faces();
    const labelListList pointPoints = mesh.pointPoints();
    const labelListList cellPoints = mesh.cellPoints();
    const labelList faceOwner = mesh.faceOwner();
    const vectorField faceNormals = mesh.Sf()/mesh.magSf();

    labelListList abaqusCellPoints(cellPoints.size(), List<label>(8, 1));

    forAll(cells, celli)
    {
        //- face0 will be our reference face
        //- face0 seems to be always owned by the cell
        //- so the face normal points out of the cell
        label refFace = cells[celli][0];

        //- insert first four abaqusCellPoints
        abaqusCellPoints[celli][0] = (faces[refFace][3] + 1);
        abaqusCellPoints[celli][1] = (faces[refFace][2] + 1);
        abaqusCellPoints[celli][2] = (faces[refFace][1] + 1);
        abaqusCellPoints[celli][3] = (faces[refFace][0] + 1);

        //- now find the opposite face in the cell
        //Info << "Finding oppFace" << endl << endl;
        labelList refFacePoints = faces[refFace];
        //- compare each faces points to the refFace, the opposite
        //- face will share points with the refFace
        label oppFace = -1;

        //- for all the cell faces
        forAll(cells[celli], facei)
        {
            bool faceHasNoCommonPoints = true;
            label globalFaceLabel = cells[celli][facei];

            //- for all the face points
            forAll(faces[globalFaceLabel], pointi)
            {
                label globalPointLabel = faces[globalFaceLabel][pointi];

                //- compare each point with all the refFace points
                forAll(faces[refFace], refFacePointi)
                {
                    label refFaceGlobalPointLabel = faces[refFace][refFacePointi];
                    if(globalPointLabel == refFaceGlobalPointLabel)
                      {
                        faceHasNoCommonPoints = false;
                      }
                }
            }
            if(faceHasNoCommonPoints)
            {
                oppFace = globalFaceLabel;
                break;
            }
        }
        if(oppFace == -1)
        {
            SeriousError << "\n\nCannot find face opposite reference face,"
                << " are you sure this is a hex mesh?"
                << endl
                << exit(FatalError);
        }


        //- Now find out the which point on oppFace is attached
        //- to each point of the refFace

       //- for all the oppFace points
       forAll(faces[oppFace], pointi)
       {
           label globalPointi = faces[oppFace][pointi];

           //- get the oppFace pointPoints
           const labelList& oppFacePPs = pointPoints[globalPointi];

           bool ppFound = false;

           //- for all the oppFace point points
           forAll(oppFacePPs, oppFacePointi)
           {
               label globalPpi = oppFacePPs[oppFacePointi];
               if(globalPpi == faces[refFace][0])
               {
                   abaqusCellPoints[celli][7] = globalPointi + 1;
                   ppFound = true;
                   break;
               }
               else if(globalPpi == faces[refFace][1])
               {
                   abaqusCellPoints[celli][6] = globalPointi + 1;
                   ppFound = true;
                   break;
               }
               else if(globalPpi == faces[refFace][2])
               {
                    abaqusCellPoints[celli][5] = globalPointi + 1;
                    ppFound = true;
                    break;
               }
               else if(globalPpi == faces[refFace][3])
               {
                   abaqusCellPoints[celli][4] = globalPointi + 1;
                   ppFound = true;
                   break;
               }
           }
           if(!ppFound)
           {
                SeriousError << "\n\nCannot find point point opposite reference face points,"
                    << " some thing strange is happening..."
                    << endl
                    << exit(FatalError);
           }
        }
    }


    //----------------------------------------------------------------------//
    //- Writing elements
    //----------------------------------------------------------------------//
    Info << "Writing Elements" << endl << endl;
    abaqusInp << "*Element, type=C3D8R" << endl;
    forAll(cells, celli)
    {
       //- print cell number
       abaqusInp << " " << (celli+1);

       //- print points on the first face
       forAll(abaqusCellPoints[celli], pointi)
       {
          abaqusInp << ", " << (abaqusCellPoints[celli][pointi]);
       }
       abaqusInp << endl;
    }


    //----------------------------------------------------------------------//
    //- Writing node sets for each patch
    //----------------------------------------------------------------------//
    Info << "Writing node sets for each boundary patch" << endl;
    forAll(mesh.boundaryMesh(), patchi)
    {
        Info << "\tWriting points on patch " << mesh.boundaryMesh()[patchi].name()
            << " to nodeSet" << mesh.boundaryMesh()[patchi].name() << endl;

        //- Node set name
        abaqusInp << "*Nset, nset=nodeSet" << mesh.boundaryMesh()[patchi].name() << endl;

        //- Node labels
        forAll(mesh.boundaryMesh()[patchi].meshPoints(), pointi)
        {
            label globalPointi = mesh.boundaryMesh()[patchi].meshPoints()[pointi];
            abaqusInp << globalPointi + 1;
            //- a comma after every label except the last
            if( pointi != (mesh.boundaryMesh()[patchi].meshPoints().size()-1) )
            {
                abaqusInp << ", ";
            }
            else
            {
                abaqusInp << endl;
            }

            //- put 10 labels on each line for tidiness
            if(pointi % 10 == 0 && pointi != 0 && pointi != (mesh.boundaryMesh()[patchi].meshPoints().size()-1))
            {
               abaqusInp << endl;
            }
        }
    }
    Info << endl;


    //----------------------------------------------------------------------//
    //- Writing element sets for each patch
    //----------------------------------------------------------------------//
    Info << "Writing element sets for each boundary patch" << endl;
    forAll(mesh.boundaryMesh(), patchi)
    {
        Info << "\tWriting cells on patch " << mesh.boundaryMesh()[patchi].name()
            << " to elementSet" << mesh.boundaryMesh()[patchi].name() << endl;

        //- Element set name
        abaqusInp << "*Elset, elset=elementSet" << mesh.boundaryMesh()[patchi].name() << endl;

        //- Element labels
        forAll(mesh.boundaryMesh()[patchi].faceCells(), celli)
        {
             label globalCelli = mesh.boundaryMesh()[patchi].faceCells()[celli];
             abaqusInp << globalCelli + 1;
             //- a comma after every label except the last
             if( celli != (mesh.boundaryMesh()[patchi].size()-1) )
             {
                 abaqusInp << ", ";
             }
             else
             {
                 abaqusInp << endl;
             }

             //- put 10 labels on each line for tidiness
             if(celli % 10 == 0 && celli != 0 && celli != (mesh.boundaryMesh()[patchi].size()-1))
             {
                 abaqusInp << endl;
             }
        }
    }
    Info << endl;


    //----------------------------------------------------------------------//
    //- Writing element set for each material
    //----------------------------------------------------------------------//
    IOobject materialsHeader
    (
        "materials",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ
    );

    // Check U exists
    if (materialsHeader.headerOk())
    {
        Info << "Reading materials field and writing each material to an element set"
            << endl;

        volScalarField materials(materialsHeader, mesh);
        const scalarField& materialsI = materials.internalField();

        //- order the material into groups and record their cell numbers
        label numberOfMaterials = round(max(materials).value()) + 1;
        labelList cellsOfEachMaterial(numberOfMaterials,0);
        forAll(materialsI, celli)
        {
            label matNum = materialsI[celli];
            forAll(cellsOfEachMaterial, mati)
            {
                if(matNum == mati)
                {
                    cellsOfEachMaterial[mati]++;
                }
            }
        }

        labelListList materialsOrdered(numberOfMaterials, List<label>(0,0));
        forAll(materialsOrdered, matOrderi)
        {
             materialsOrdered[matOrderi].setSize(cellsOfEachMaterial[matOrderi], 0);
             label mc = 0;
             forAll(materialsI, celli)
             {
                 label matNum = materialsI[celli];
                 if(matNum == matOrderi)
                 {
                     materialsOrdered[matOrderi][mc] = celli;
                     mc++;
                }
             }

             Info << "\tWriting " << cellsOfEachMaterial[matOrderi] << " cells of material "
                 << matOrderi << " to elementSetMaterial" << matOrderi << endl;
             abaqusInp << "*Elset, elset=elementSetMaterial" << matOrderi << endl;
             forAll(materialsOrdered[matOrderi], celli)
             {
                 abaqusInp << (materialsOrdered[matOrderi][celli]+1);
                 if(celli != materialsOrdered[matOrderi].size()-1)
                 {
                     abaqusInp << "," << endl;
                 }
                 else
                 {
                     abaqusInp << endl;
                }
           }
        }
    }
    else
    {
        Info << "No materials file present so no material element sets will be written"
            << endl << endl;
    }
    Info << endl;


    //----------------------------------------------------------------------//
    //- Writing surfaces for each patch
    //----------------------------------------------------------------------//
    Info << "Writing surfaces for each boundary patch" << endl;
    forAll(mesh.boundaryMesh(), patchi)
    {
        Info << "\tWriting cells on patch " << mesh.boundaryMesh()[patchi].name()
            << " to surface" << mesh.boundaryMesh()[patchi].name() << endl;

        //- Surface name
        abaqusInp << "*Surface, type=ELEMENT, name=surface" << mesh.boundaryMesh()[patchi].name() << endl;

        //- the face is denoted by either s1, s2, s3, s4, s5 or s6
        //- so I have to find the exterior face of the element and
        //- the corresponding s number
        //- face convention
        //- face 1 is nodes 1 2 3 4
        //- face 2 is nodes 5 8 7 6
        //- face 3 is nodes 1 5 6 2
        //- face 4 is nodes 2 6 7 3
        //- face 5 is nodes 3 7 8 4
        //- face 6 is nodes 4 8 5 1
        //- So i will compare points on the face and
        labelListList abaqusFaceConvention(6, List<label>(4,-1));
        abaqusFaceConvention[0][0] = 1;
        abaqusFaceConvention[0][1] = 2;
        abaqusFaceConvention[0][2] = 3;
        abaqusFaceConvention[0][3] = 4;
        abaqusFaceConvention[1][0] = 5;
        abaqusFaceConvention[1][1] = 8;
        abaqusFaceConvention[1][2] = 7;
        abaqusFaceConvention[1][3] = 6;
        abaqusFaceConvention[2][0] = 1;
        abaqusFaceConvention[2][1] = 5;
        abaqusFaceConvention[2][2] = 6;
        abaqusFaceConvention[2][3] = 2;
        abaqusFaceConvention[3][0] = 2;
        abaqusFaceConvention[3][1] = 6;
        abaqusFaceConvention[3][2] = 7;
        abaqusFaceConvention[3][3] = 3;
        abaqusFaceConvention[4][0] = 3;
        abaqusFaceConvention[4][1] = 7;
        abaqusFaceConvention[4][2] = 8;
        abaqusFaceConvention[4][3] = 4;
        abaqusFaceConvention[5][0] = 4;
        abaqusFaceConvention[5][1] = 8;
        abaqusFaceConvention[5][2] = 5;
        abaqusFaceConvention[5][3] = 1;

        //- Element labels and the free face label
        forAll(mesh.boundaryMesh()[patchi], facei)
        {
            label globalCelli = mesh.boundaryMesh()[patchi].faceCells()[facei];
            abaqusInp << (globalCelli + 1) << ", ";

            //- compare this face's points with the abaqusCellPoints and see which points
            //- are used, then compared the points used to the faces above to see which
            //- face it is
            const labelList& facePoints = mesh.boundaryMesh()[patchi][facei];
            const labelList& thisCellAbaqusPoints = abaqusCellPoints[globalCelli];
            labelList thisFaceAbaqusPoints(4,-1);
            label pointsFound = 0;
            forAll(thisCellAbaqusPoints, pointi)
            {
                label cellPointi = thisCellAbaqusPoints[pointi];
                if(cellPointi == facePoints[0]+1)
                {
                    pointsFound++;
                    thisFaceAbaqusPoints[0] = pointi+1;
                }
                else if(cellPointi == facePoints[1]+1)
                {
                    pointsFound++;
                    thisFaceAbaqusPoints[1] = pointi+1;
                }
                else if(cellPointi == facePoints[2]+1)
                {
                    pointsFound++;
                    thisFaceAbaqusPoints[2] = pointi+1;
                }
                else if(cellPointi == facePoints[3]+1)
                {
                    pointsFound++;
                    thisFaceAbaqusPoints[3] = pointi+1;
                }
            }
            if(pointsFound != 4)
            {
                SeriousError << "Something is wrong while determing the surfaces from the patches"
                    << endl << exit(FatalError);
            }

            //- now thisFaceAbaqusPoints holds the points on the face so
            //- compare this face to face[1-6] to see which face it is.
            //- The points won't neccesarily be in the same order as faces[1-6]
            //- but will contain the same labels
            forAll(abaqusFaceConvention, abqfacei)
            {
                int pointsOnFace = 0;
                forAll(abaqusFaceConvention[abqfacei], pointi)
                {
                    forAll(thisFaceAbaqusPoints, pi)
                    {
                        if(thisFaceAbaqusPoints[pi] == abaqusFaceConvention[abqfacei][pointi])
                        {
                            pointsOnFace++;
                        }
                    }
                }
                if(pointsOnFace == 4)
                {
                    //Info << "found face!" << endl;
                    abaqusInp << " s" << (abqfacei + 1) << endl;
                }
            }
        }
    }
    Info << endl;


    //----------------------------------------------------------------------//
    //- end part
    //----------------------------------------------------------------------//
    abaqusInp << "*End Part" << endl;

    Info << "Finished writing abaqusMesh.inp" << endl << endl;

    return(0);
}


// ************************************************************************* //