520 lines
18 KiB
C++
520 lines
18 KiB
C++
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / 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 <http://www.gnu.org/licenses/>.
|
|
|
|
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);
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|