/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | foam-extend: Open Source CFD \\ / O peration | Version: 4.1 \\ / A nd | Web: http://www.foam-extend.org \\/ M anipulation | For copyright notice see file Copyright ------------------------------------------------------------------------------- License This file is part of foam-extend. foam-extend 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 3 of the License, or (at your option) any later version. foam-extend 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 foam-extend. If not, see . Application foamMeshToAbaqus Description Converts an FOAM 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(), labelList(8, label(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, labelList(0, label(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, labelList(4, label(-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, label(-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); } // ************************************************************************* //