diff --git a/libsrc/interface/writeabaqus.cpp b/libsrc/interface/writeabaqus.cpp index 8a0557c9..8799ef0f 100644 --- a/libsrc/interface/writeabaqus.cpp +++ b/libsrc/interface/writeabaqus.cpp @@ -15,15 +15,103 @@ namespace netgen { +using std::vector; +struct AbaqusElementType +{ + const char * name; + const vector permutation; + AbaqusElementType(const char * name, const vector & permutation) + : name(name), permutation(permutation) + {} +}; + +static inline const AbaqusElementType & GetAbaqusType(int dim, int num_nodes) +{ + // maps num_nodes to AbaqusElementType for each dimension + typedef std::map AbaqusElementTypes; + static const std::map abaqus_eltypes[3] = + { + // 1D + AbaqusElementTypes{ + {2, AbaqusElementType{"T2D2", vector{0,1}}}, + }, + // 2D + AbaqusElementTypes{ + {3, AbaqusElementType{"CPS3", vector{0,1,2}}}, + }, + // 3D + AbaqusElementTypes{ + {4, AbaqusElementType{"C3D4", vector{0,1,3,2}}}, + {10, AbaqusElementType{"C3D10", vector{0,1,3,2,4,8,6,5,7,9}}}, + } + }; + + const auto & eltypes = abaqus_eltypes[dim-1]; + if (eltypes.count(num_nodes) > 0) + return eltypes.at(num_nodes); + else + throw Exception("unsupported " + ToString(dim)+"d Element type with " + ToString(num_nodes) + " nodes"); +} + +static void WritePoints ( const Mesh & mesh, ostream & out ) +{ + out << "*Node" << endl; + for(auto pi : mesh.Points().Range() ) + { + out << pi+1-PointIndex::BASE << ", "; + auto p = mesh[pi]; + out << p[0] << ", " << p[1] << ", " << p[2] << '\n'; + } +} + +template +static void WriteElement(ostream & out, const Mesh& mesh, ElIndex ei, const vector & permutation, int & el_counter) +{ + el_counter++; + auto el = mesh[ei]; + out << el_counter; + for(auto i : Range(el.PNums())) + out << ", " << el[permutation[i]]+1-PointIndex::BASE; + out << '\n'; +} + +template +static void WriteElements ( ostream & out, const Mesh & mesh, int dim, const Elements & el_range, int & el_counter) +{ + // map index, num_nodes to elements + std::map, Array> elset_map; + + for(auto ei : el_range) + { + const auto & el = mesh[ei]; + int index = 0; + if constexpr(std::is_same_v) + index = el.edgenr; + else + index = el.GetIndex(); + elset_map[{index, el.GetNP()}].Append(ei); + } + + for(auto & [key, elems] : elset_map) + { + auto [index, num_nodes] = key; + auto name = mesh.GetRegionName(elems[0]); + if (name == "") name = "default"; + PrintMessage (5, index, ": ", name); + const auto & eltype = GetAbaqusType(dim, num_nodes) ; + out << "*Element, type=" << eltype.name << ", ELSET=" << name << endl; + for(auto ei : elems) + WriteElement(out, mesh, ei, eltype.permutation, el_counter); + } +} void WriteAbaqusFormat (const Mesh & mesh, const filesystem::path & filename) { - - cout << "\nWrite Abaqus Volume Mesh" << endl; + PrintMessage (1, "Write Abaqus Mesh"); ofstream outfile (filename); @@ -32,96 +120,17 @@ void WriteAbaqusFormat (const Mesh & mesh, outfile.precision(8); - outfile << "*Node" << endl; - - int np = mesh.GetNP(); - int ne = mesh.GetNE(); - 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"; - } - - 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;} - } + int element_counter = 0; + WritePoints(mesh, outfile); + if(mesh.GetDimension() < 3) + WriteElements(outfile, mesh, 1, mesh.LineSegments().Range(), element_counter); + WriteElements(outfile, mesh, 2, mesh.SurfaceElements().Range(), element_counter); + WriteElements(outfile, mesh, 3, mesh.VolumeElements().Range(), element_counter); + // Write identifications (untested!) if (mesh.GetIdentifications().GetMaxNr()) { + const auto np = mesh.GetNP(); // periodic identification, implementation for // Helmut J. Boehm, TU Vienna @@ -138,27 +147,27 @@ void WriteAbaqusFormat (const Mesh & mesh, NgArray pairs; NgBitArray master(np), help(np); master.Set(); - for (i = 1; i <= 3; i++) + for (int i = 1; i <= 3; i++) { mesh.GetIdentifications().GetPairs (i, pairs); help.Clear(); - for (j = 1; j <= pairs.Size(); j++) + for (int j = 1; j <= pairs.Size(); j++) { help.Set (pairs.Get(j).I1()); } master.And (help); } - for (i = 1; i <= np; i++) + for (int i = 1; i <= np; i++) if (master.Test(i)) masternode = i; cout << "masternode = " << masternode << " = " << mesh.Point(masternode) << endl; NgArray minions(3); - for (i = 1; i <= 3; i++) + for (int i = 1; i <= 3; i++) { mesh.GetIdentifications().GetPairs (i, pairs); - for (j = 1; j <= pairs.Size(); j++) + for (int j = 1; j <= pairs.Size(); j++) { if (pairs.Get(j).I1() == masternode) minions.Elem(i) = pairs.Get(j).I2(); @@ -179,12 +188,12 @@ void WriteAbaqusFormat (const Mesh & mesh, << "**POINT_fixed\n" << "**\n" << "*BOUNDARY, OP=NEW\n"; - for (j = 1; j <= 3; j++) + for (int j = 1; j <= 3; j++) outfile << masternode << ", " << j << ",, 0.\n"; outfile << "**\n" << "*BOUNDARY, OP=NEW\n"; - for (j = 1; j <= 3; j++) + for (int j = 1; j <= 3; j++) { Vec3d v(mesh.Point(masternode), mesh.Point(minions.Get(j))); double vlen = v.Length(); @@ -203,18 +212,18 @@ void WriteAbaqusFormat (const Mesh & mesh, NgBitArray eliminated(np); eliminated.Clear(); - for (i = 1; i <= mesh.GetIdentifications().GetMaxNr(); i++) + for (int i = 1; i <= mesh.GetIdentifications().GetMaxNr(); i++) { mesh.GetIdentifications().GetPairs (i, pairs); if (!pairs.Size()) continue; - for (j = 1; j <= pairs.Size(); j++) + for (int 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++) + for (int k = 1; k <= 3; k++) { mpc << "4" << "\n"; mpc << pairs.Get(j).I2() << "," << k << ", -1.0, "; @@ -227,7 +236,7 @@ void WriteAbaqusFormat (const Mesh & mesh, } - cout << "done" << endl; + PrintMessage(1, "done"); } static RegisterUserFormat reg_abaqus ("Abaqus Format", {".mesh"}, nullopt, WriteAbaqusFormat);