mirror of
https://github.com/NGSolve/netgen.git
synced 2025-05-19 00:30:48 +05:00
348 lines
9.9 KiB
C++
348 lines
9.9 KiB
C++
//
|
|
// Write Abaqus file
|
|
//
|
|
//
|
|
|
|
#include <mystdlib.h>
|
|
|
|
#include <myadt.hpp>
|
|
#include <linalg.hpp>
|
|
#include <csg.hpp>
|
|
#include <meshing.hpp>
|
|
|
|
namespace netgen
|
|
{
|
|
#include "writeuser.hpp"
|
|
|
|
|
|
void WriteAbaqusVolumeElements(const Mesh & mesh,
|
|
ofstream &outfile) {
|
|
int i, j, k;
|
|
int ne = mesh.GetNE();
|
|
int elemcnt = 0; //element counter
|
|
int finished = 0;
|
|
int indcnt = 1; //index counter
|
|
|
|
while (!finished)
|
|
{
|
|
int actcnt = 0;
|
|
const Element & el1 = mesh.VolumeElement(1);
|
|
int non = el1.GetNP();
|
|
if (non == 4)
|
|
{
|
|
outfile << "*Element, type=C3D4, ELSET=PART" << indcnt << endl;
|
|
}
|
|
else if (non == 10)
|
|
{
|
|
outfile << "*Element, type=C3D10, ELSET=PART" << indcnt << endl;
|
|
}
|
|
else
|
|
{
|
|
cout << "unsupported Element type!!!" << endl;
|
|
}
|
|
|
|
for (i = 1; i <= ne; i++)
|
|
{
|
|
const Element & el = mesh.VolumeElement(i);
|
|
|
|
if (el.GetIndex() == indcnt)
|
|
{
|
|
actcnt++;
|
|
if (el.GetNP() != non)
|
|
{
|
|
cout << "different element-types in a subdomain are not possible!!!" << endl;
|
|
continue;
|
|
}
|
|
|
|
elemcnt++;
|
|
outfile << elemcnt << ", ";
|
|
if (non == 4)
|
|
{
|
|
outfile << el.PNum(1) << ", ";
|
|
outfile << el.PNum(2) << ", ";
|
|
outfile << el.PNum(4) << ", ";
|
|
outfile << el.PNum(3) << "\n";
|
|
}
|
|
else if (non == 10)
|
|
{
|
|
outfile << el.PNum(1) << ", ";
|
|
outfile << el.PNum(2) << ", ";
|
|
outfile << el.PNum(4) << ", ";
|
|
outfile << el.PNum(3) << ", ";
|
|
outfile << el.PNum(5) << ", ";
|
|
outfile << el.PNum(9) << ", ";
|
|
outfile << el.PNum(7) << ", " << "\n";
|
|
outfile << el.PNum(6) << ", ";
|
|
outfile << el.PNum(8) << ", ";
|
|
outfile << el.PNum(10) << "\n";
|
|
}
|
|
else
|
|
{
|
|
cout << "unsupported Element type!!!" << endl;
|
|
for (j = 1; j <= el.GetNP(); j++)
|
|
{
|
|
outfile << el.PNum(j);
|
|
if (j != el.GetNP()) outfile << ", ";
|
|
}
|
|
outfile << "\n";
|
|
}
|
|
}
|
|
}
|
|
indcnt++;
|
|
if (elemcnt == ne) {finished = 1; cout << "all elements found by Index!" << endl;}
|
|
if (actcnt == 0) {finished = 1;}
|
|
}
|
|
}
|
|
|
|
void WriteAbaqusSurfaceElements(const Mesh & mesh,
|
|
ofstream& outfile) {
|
|
int i, j, k;
|
|
int ne = mesh.GetNSE();
|
|
int elemcnt = 0; //element counter
|
|
int finished = 0;
|
|
int indcnt = 1; //index counter
|
|
|
|
while (!finished)
|
|
{
|
|
int actcnt = 0;
|
|
const Element2d & el1 = mesh.SurfaceElement(1);
|
|
int non = el1.GetNP();
|
|
if (non == 3)
|
|
{
|
|
outfile << "*Element, type=S3, ELSET=PART" << indcnt << endl;
|
|
}
|
|
else if (non == 6)
|
|
{
|
|
outfile << "*Element, type=S6, ELSET=PART" << indcnt << endl;
|
|
}
|
|
else
|
|
{
|
|
cout << "unsupported Element type!!!" << endl;
|
|
}
|
|
|
|
for (i = 1; i <= ne; i++)
|
|
{
|
|
const Element2d & el = mesh.SurfaceElement(i);
|
|
|
|
if (el.GetIndex() == indcnt)
|
|
{
|
|
actcnt++;
|
|
if (el.GetNP() != non)
|
|
{
|
|
cout << "different element-types in a subdomain are not possible!!!" << endl;
|
|
continue;
|
|
}
|
|
|
|
elemcnt++;
|
|
outfile << elemcnt << ", ";
|
|
if (non == 3)
|
|
{
|
|
outfile << el.PNum(1) << ", ";
|
|
outfile << el.PNum(2) << ", ";
|
|
outfile << el.PNum(3) << "\n";
|
|
}
|
|
else if (non == 6)
|
|
{
|
|
outfile << el.PNum(1) << ", ";
|
|
outfile << el.PNum(2) << ", ";
|
|
outfile << el.PNum(3) << ", ";
|
|
outfile << el.PNum(4) << ", ";
|
|
outfile << el.PNum(5) << ", ";
|
|
outfile << el.PNum(6) << "\n";
|
|
}
|
|
else
|
|
{
|
|
cout << "unsupported Element type!!!" << endl;
|
|
for (j = 1; j <= el.GetNP(); j++)
|
|
{
|
|
outfile << el.PNum(j);
|
|
if (j != el.GetNP()) outfile << ", ";
|
|
}
|
|
outfile << "\n";
|
|
}
|
|
}
|
|
}
|
|
indcnt++;
|
|
if (elemcnt == ne) {finished = 1; cout << "all elements found by Index!" << endl;}
|
|
if (actcnt == 0) {finished = 1;}
|
|
}
|
|
}
|
|
|
|
|
|
void WriteAbaqusFormat (const Mesh & mesh,
|
|
const string & filename,
|
|
bool forcesurface,
|
|
bool forcevolume)
|
|
|
|
{
|
|
|
|
int nve = mesh.GetNE();
|
|
int nse = mesh.GetNSE();
|
|
|
|
bool writevolume = forcevolume || nve>0;
|
|
bool writesurface = forcesurface || (nve == 0 && nse>0);
|
|
|
|
if ( writevolume && writesurface)
|
|
{
|
|
cout << "\nWrite Abaqus Volume and Surface Mesh" << endl;
|
|
}
|
|
else if ( writevolume )
|
|
{
|
|
cout << "\nWrite Abaqus Volume Mesh" << endl;
|
|
}
|
|
else if ( writesurface )
|
|
{
|
|
cout << "\nWrite Abaqus Surface Mesh" << endl;
|
|
}
|
|
else
|
|
{
|
|
cout << "\nWrite Abaqus Mesh Nodes" << endl;
|
|
}
|
|
|
|
ofstream outfile (filename.c_str());
|
|
|
|
outfile << "*Heading" << endl;
|
|
outfile << " " << filename << endl;
|
|
|
|
outfile.precision(8);
|
|
|
|
outfile << "*Node" << endl;
|
|
|
|
int np = mesh.GetNP();
|
|
int i, j, k;
|
|
|
|
for (i = 1; i <= np; i++)
|
|
{
|
|
outfile << i << ", ";
|
|
outfile << mesh.Point(i)(0) << ", ";
|
|
outfile << mesh.Point(i)(1) << ", ";
|
|
outfile << mesh.Point(i)(2) << "\n";
|
|
}
|
|
|
|
if ( forcevolume || nve>0 )
|
|
{
|
|
WriteAbaqusVolumeElements(mesh, outfile);
|
|
}
|
|
if ( forcesurface || (nse>0 && nve == 0) )
|
|
{
|
|
WriteAbaqusSurfaceElements(mesh, outfile);
|
|
}
|
|
|
|
|
|
if (mesh.GetIdentifications().GetMaxNr())
|
|
{
|
|
// periodic identification, implementation for
|
|
// Helmut J. Boehm, TU Vienna
|
|
|
|
char cfilename[255];
|
|
strcpy (cfilename, filename.c_str());
|
|
|
|
char mpcfilename[255];
|
|
strcpy (mpcfilename, cfilename);
|
|
size_t len = strlen (cfilename);
|
|
if (len >= 4 && (strcmp (mpcfilename+len-4, ".inp") == 0))
|
|
strcpy (mpcfilename+len-4, ".mpc");
|
|
else
|
|
strcat (mpcfilename, ".mpc");
|
|
|
|
ofstream mpc (mpcfilename);
|
|
|
|
int masternode(0);
|
|
|
|
NgArray<INDEX_2> pairs;
|
|
NgBitArray master(np), help(np);
|
|
master.Set();
|
|
for (i = 1; i <= 3; i++)
|
|
{
|
|
mesh.GetIdentifications().GetPairs (i, pairs);
|
|
help.Clear();
|
|
for (j = 1; j <= pairs.Size(); j++)
|
|
{
|
|
help.Set (pairs.Get(j).I1());
|
|
}
|
|
master.And (help);
|
|
}
|
|
for (i = 1; i <= np; i++)
|
|
if (master.Test(i))
|
|
masternode = i;
|
|
|
|
cout << "masternode = " << masternode << " = "
|
|
<< mesh.Point(masternode) << endl;
|
|
NgArray<int> slaves(3);
|
|
for (i = 1; i <= 3; i++)
|
|
{
|
|
mesh.GetIdentifications().GetPairs (i, pairs);
|
|
for (j = 1; j <= pairs.Size(); j++)
|
|
{
|
|
if (pairs.Get(j).I1() == masternode)
|
|
slaves.Elem(i) = pairs.Get(j).I2();
|
|
}
|
|
cout << "slave(" << i << ") = " << slaves.Get(i)
|
|
<< " = " << mesh.Point(slaves.Get(i)) << endl;
|
|
}
|
|
|
|
|
|
outfile << "**\n"
|
|
<< "*NSET,NSET=CTENODS\n"
|
|
<< slaves.Get(1) << ", "
|
|
<< slaves.Get(2) << ", "
|
|
<< slaves.Get(3) << endl;
|
|
|
|
|
|
outfile << "**\n"
|
|
<< "**POINT_fixed\n"
|
|
<< "**\n"
|
|
<< "*BOUNDARY, OP=NEW\n";
|
|
for (j = 1; j <= 3; j++)
|
|
outfile << masternode << ", " << j << ",, 0.\n";
|
|
|
|
outfile << "**\n"
|
|
<< "*BOUNDARY, OP=NEW\n";
|
|
for (j = 1; j <= 3; j++)
|
|
{
|
|
Vec3d v(mesh.Point(masternode), mesh.Point(slaves.Get(j)));
|
|
double vlen = v.Length();
|
|
int dir = 0;
|
|
if (fabs (v.X()) > 0.9 * vlen) dir = 2;
|
|
if (fabs (v.Y()) > 0.9 * vlen) dir = 3;
|
|
if (fabs (v.Z()) > 0.9 * vlen) dir = 1;
|
|
if (!dir)
|
|
cout << "ERROR: Problem with rigid body constraints" << endl;
|
|
outfile << slaves.Get(j) << ", " << dir << ",, 0.\n";
|
|
}
|
|
|
|
outfile << "**\n"
|
|
<< "*EQUATION, INPUT=" << mpcfilename << endl;
|
|
|
|
|
|
NgBitArray eliminated(np);
|
|
eliminated.Clear();
|
|
for (i = 1; i <= mesh.GetIdentifications().GetMaxNr(); i++)
|
|
{
|
|
mesh.GetIdentifications().GetPairs (i, pairs);
|
|
if (!pairs.Size())
|
|
continue;
|
|
|
|
for (j = 1; j <= pairs.Size(); j++)
|
|
if (pairs.Get(j).I1() != masternode &&
|
|
!eliminated.Test(pairs.Get(j).I2()))
|
|
{
|
|
eliminated.Set (pairs.Get(j).I2());
|
|
for (k = 1; k <= 3; k++)
|
|
{
|
|
mpc << "4" << "\n";
|
|
mpc << pairs.Get(j).I2() << "," << k << ", -1.0, ";
|
|
mpc << pairs.Get(j).I1() << "," << k << ", 1.0, ";
|
|
mpc << slaves.Get(i) << "," << k << ", 1.0, ";
|
|
mpc << masternode << "," << k << ", -1.0 \n";
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
cout << "done" << endl;
|
|
}
|
|
|
|
}
|